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Introduction: 

Recently  optical  diagnostics  based  on  diffuse  near-infrared  (NIR)  light  has  been  employed  in 
breast  cancer  detection  [1-4].  Optical  spectroscopy  and  imaging  of  body  structure  and  function  is 
made  possible  by  a  spectral  window  that  exists  within  tissues  in  the  700  -  900  nm  region  (near- 
infrared  (NIR)),  in  which  photon  transport  is  dominated  by  scattering  rather  than  absorption. 
Thus,  to  a  very  good  approximation,  NIR  photons  diffuse  through  relatively  thick  tissues. 
Compared  with  other  imaging  methods,  the  optical  method  has  numerous  advantages.  NIR  light 
has  high  specificity  because  certain  characteristics  differ  in  tumors  and  normal  breast  tissues:  1) 
The  amount  of  blood  needed  to  serve  the  tumor’s  metabolic  needs  is  increased  over  that  of  normal 
background  tissue;  2)  Hemoglobin  desaturation  in  tumors  is  increased  due  to  the  high  oxygen 
demand  of  cancers  [5];  3)  Light  scattering  of  the  tumor  is  enhanced  due  to  the  increased 

mitochondria  population  with  respect  to  the  background  of  normal  cells  [6];  4)  The  angiogenesis 
phenomenon  [7]  allows  tumor  identification  by  observing  the  delivery  of  contrast  agents,  such  as 
the  well-characterized  NIR  absorber  indocyanine  green  (ICG),  through  the  permeable  vascular 
beds  to  the  extravascular  space  surrounding  the  tumor.  In  addition,  the  NIR  device  requires 
relatively  low  manufacturing  and  operating  costs,  NIR  light  is  non-invasive  to  the  human  body, 
and  the  device  is  compact  and  easily  portable.  As  a  potential  diagnostic,  however,  diffused  NIR 
light  imaging  and  spectroscopy  suffer  from  low  spatial  resolution  due  to  the  diffusive  nature  of 
photon  density  waves.  Current  instruments  can  distinguish  simple  structures  of  approximately  1 
cm  in  size,  but  sharp  edges  are  typically  blurred  by  a  few  millimeters  [2],  In  addition,  imaging 
reconstruction  is  based  on  inverse  scattering  approaches  [8-10]  and  the  inverse  problem  is  in 
genera]  underdetermined  and  ill-posed.  Thus  while  the  optical  method  offers  new  routes  to 
tumor  specificity,  the  relatively  low  resolution  and  difficulties  in  image  reconstruction  make  it 
difficult  to  take  full  advantage  of  this  contrast. 

The  concept  of  pulse-echo  ultrasound  imaging  is  well  known.  Ultrasound  has  many  advantages: 
its  manufacturing  and  operating  costs  are  low;  it  is  non-invasive  to  the  human  body  and  it 
provides  high-resolution  real-time  images.  Current  state-of-the-art  breast  scanners  can  detect 
small  lesions  of  several  millimeters  in  size  [13].  Ultrasound  is  frequently  used  in  conjunction 
with  mammography  to  differentiate  simple  cysts  from  solid  lesions.  When  the  criteria  for  a 
simple  cyst  are  strictly  adhered  to,  the  accuracy  of  ultrasound  is  96%-100%.  However,  the 
ultrasound  appearance  of  benign  and  malignant  lesions  has  considerable  overlapping  features  [11- 
12],  which  has  prompted  many  radiologists  to  recommend  biopsies  on  most  solid  nodules.  This 
results  in  a  large  number  of  biopsies  yielding  normal  or  benign  breast  tissues. 

We  have  developed  a  novel  hybrid  imaging  technique,  which  combines  the  high  contrasts  of  NIR 
imaging  in  distinguishing  between  normal,  benign  and  malignant  tissues  with  the  high  spatial 
resolution  inherent  in  ultrasound  imaging.  This  hybrid  imaging  method  provides  a  new  way  to 
utilize  optical  contrast  and  ultrasound  imaging  capability  [14-17]. 

In  this  DOD  ARMY  sponsored  study,  we  proposed  to  achieve  the  flowing  objectives 

1)  to  refine  our  existing  NIR  optical  imaging  system  hardware  toward  high  signal-to-noise  ratio 
and  fast  data  acquisition  (task  1); 

2)  to  implement  to  imaging  software  for  clinical  studies  (task  1) 

3)  to  optimize  the  combined  probe  design  through  simulation  and  phantom  experiments  (task  1); 

4)  to  validate  the  combined  imaging  in  cancer  detection  and  diagnosis  through  clinical  studies 
(task  2). 


zhu 


Page  3 


9/16/2001 


We  have  successfully  completed  task  1  proposed  in  the  application  and  has  started  the  task  2  by 
recruiting  breast  cancer  patients.  We  expect  to  start  clinical  studies  in  a  month. 

Body 

Under  the  supports  of  DOD  ARMY  we  have  completed  the  combined  probe  by  co-axially 
deploying  optical  and  ultrasound  probes  simultaneously.  This  combined  probe  can  partially 
overcome  the  lesion  mapping  uncertainty  problem  encountered  in  our  preliminary  clinical 
studies  where  optical  and  ultrasound  probes  were  used  separately  to  scan  deformable  breasts 
[14].  The  optical  probe  consists  of  12  dual  wavelength  sources  and  8  optical  detectors  [17],  and 
the  sources  and  detectors  are  coupled  through  optical  fibers  to  a  hand-held  probe  shown  in 
Figure  1(a).  A  commercial  ultrasound  probe  (7  MHz  linear  array)  is  simultaneously  deployed  in 
the  middle  of  the  combined  probe.  With  this  combined  probe,  we  can  demonstrate  the  advantage 
of  dual-modality  diagnosis.  The  picture  of  our  NIR  scanner  is  shown  in  Figure  1(b).  All  the 
imaging  algorithms  are  implemented  in  the  software  and  the  system  will  be  used  to  perform  pilot 
clinical  studies. 


Figure  1.  (a).  Our  combined  probe  with  NIR  sensors  deployed  at  the  periphery  and  a  1-D  commercial  ultrasound 
array  located  in  the  middle,  (b).  Commercial  ultrasound  scanner  (left),  our  NIR  system  (right)  and  the  combined 
probe  (middle). 

In  addition  to  complete  the  above  proposed  tasks,  we  have  also  done  three  dimensional 
simultaneous  ultrasound  and  NIR  imaging  and  co-registration  in  phantom  tests.  In  the  on-going 
clinical  studies,  conventional  ultrasound  images  are  obtained  in  x-z  planes  and  NIR  images  are  in 
x-y  planes.  Therefore,  the  ultrasound  and  NIR  imaging  planes  are  orthogonal  to  each  other. 
Early  detection  of  breast  cancers  requires  simultaneous  ultrasound  and  NIR  imaging  at  the  same 
imaging  plane,  referred  to  as  co-registration  in  the  proposal.  Co-registration  enables  the  use 
of  lesion  morphology  provided  by  high-resolution  ultrasound  to  improve  the  lesion  optical 
property  estimate.  This  task  can  be  accomplished  by  simultaneously  deploying  a  2-D  ultrasound 
array,  which  is  capable  of  providing  3-D  high-resolution  volumetric  images,  and  the  multiple 
source  and  detector  NIR  array.  We  have  built  a  simple  2-D  ultrasound  array  consisting  of  64 
transducers,  which  were  deployed  simultaneously  with  the  NIR  source  and  detector  fibers  on  the 
same  probe,  to  demonstrate  the  feasibility  of  co-registration  [17]. 

In  Refs  [16-17],  we  have  demonstrated  that  with  the  a  priori  knowledge  of  lesion  location  and 
shape  information  provided  by  co-registered  ultrasound,  NIR  imaging  reconstruction  can  be 
localized  within  specified  spatial  and  temporal  regions.  As  a  result,  the  reconstruction  is  over- 
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determined  because  the  total  number  of  unknown  optical  properties  is  reduced  significantly.  In 
addition,  the  reconstruction  is  less  sensitive  to  noise  because  the  convergence  can  be  achieved 
within  a  small  number  of  iterations.  We  have  conducted  a  series  of  experiments  to  assess  the 
improvement  on  reconstructed  optical  absorption  coefficients  with  the  a  priori  target  depth  [16] 
and  spatial  location  information  [17].  For  the  experiments  reported  in  Ref  [16],  we  have  shown 
that  with  the  a  priori  target  depth  information  provided  by  co-registered  ultrasound,  the  accuracy 
of  the  reconstructed  absorption  coefficient  has  been  improved  by  15%  and  30%  on  average,  and 
the  resolution  measured  at  Full  Width  at  Half  Maximum  has  been  improved  by  24%  and  41%  on 
average  for  high  and  low  contrast  cases,  respectively.  The  speed  of  reconstruction  has  been 
improved  by  10  times  on  average.  For  the  experiments  reported  in  Ref  [17],  we  have 
demonstrated  that,  with  the  guidance  of  a  priori  target  temporal  and  spatial  distributions,  the 
iterative  inversion  algorithm  converges  very  fast  and  only  one  iteration  is  needed  to  obtain 
accurate  optical  absorption  coefficient.  This  result  is  very  encouraging  because  there  is  no 
known  robust  stopping  criterion  in  the  literature  to  terminate  iterative  inversion  algorithms.  With 
the  a  priori  knowledge  from  co-registered  ultrasound,  no  further  iterations  are  needed. 

Our  preliminary  results  of  using  co-registered  ultrasound  to  improve  the  accuracy  of  estimated 
lesion  optical  properties  are  very  encouraging.  However,  the  image  quality  of  our  simple  2-D 
ultrasound  array  is  far  from  commercial  standard.  We  have  upgraded  our  simple  2-D  ultrasound 
array  to  commercial  standards  by  incorporating  a  state-of-the-art  1.75D  ultrasound  array 
purchased  from  Tetrad  Inc  (see  Fig.2a)  and  constructing  the  electronic  data  acquisition  system 
(Fig.2b).  The  array  consists  of  1280  transducer  elements  and  the  ultrasound  beam  can  be 
scanned  in  both  x  and  y  spatial  directions.  Preliminary  results  obtained  from  our  high-resolution 
hybrid  imaging  system  are  given  in  Figure  3.  Figure  3(a)  is  the  ultrasound  x-z  image  of  a  small 
emulated  lesion  (6  mm  in  diameter)  of  low  optical  contrast  (absorption  coefficient  jUu  =  0. 1  cm' 

')  located  2.5  cm  deep  into  the  optical  scattering  medium  ( jus  =6 cm'1).  Figure  3(b)  is  the 
ultrasound  C-scan  or  spatial  image  of  the  target  co-registered  with  NIR  images  shown  in  Fig.3 
(c)  and  (d).  Without  ultrasound  image  guidance  shown  in  (b),  the  pa  image  shown  in  Fig.  3(c) 
does  not  converge  and  two  targets  appear  at  wrong  locations.  With  the  ultrasound  guidance 
shown  in  (b),  the  target  shown  in  Fig.3  (d)  appears  at  the  correct  location  and  the  estimated  pu  is 

86%  of  the  true  value.  These  preliminary  results  are  very  encouraging  and  suggest  that  the 
hybrid  imaging  method  has  a  great  potential  to  improve  early  detection  by  combining  high 
spatial  resolution  of  ultrasound  and  high  contrast  of  optical  imaging. 
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Figure  3.  Ultrasound  and  optical  images  of  a  small  (6mm  in  diameter)  emulated  tumor  of  low  optical 
contrast  (//a  =0.1  cm'1)  imbedded  in  the  optical  scattering  medium  (reduced  scattering  coefficient  is 
//’  =  6  cm'1),  (a)  Ultrasound  image  obtained  in  x-z  scan,  where  x  is  the  lateral  direction  and  z  is  the 

propagation  direction).  The  target  is  slightly  echogenic  and  clearly  seen,  (b)  C-scan  or  spatial  x-y 
image  of  the  ultrasound  co-registered  with  NIR  images  shown  in  (c)  and  (d).  The  horizontal  size  of 
the  image  is  3  cm  and  the  vertical  size  is  1 .7  cm.  The  size  is  scaled  to  the  same  as  in  (c)  and  (d).  (c) 
Reconstructed  optical  absorption  image  jua  without  ultrasound  guidance.  The  small  target  should 

appear  at  the  center  location  of  the  image,  however,  two  targets  appear  at  wrong  locations,  (d)  With 
the  guidance  of  co-registered  ultrasound  shown  in  (b),  the  reconstructed  target  appears  at  correct 

A  1  .  . 

location  and  the  reconstructed  value  is  fia  =0.086  cm'  which  is  86%  of  the  true  value. 


We  have  evaluated  our  new  imaging  method  with  an  excised  tumor  using  the  system  shown  in 
Fig.  lb.  The  tumor  (9L  model)  was  grown  in  a  male  Fisher  rat  and  was  excised  immediately  after 
the  rat  died  from  the  anesthesia.  The  weight  of  the  tumor  was  3.6  grams  and  the  volume  was 
about  3.6  ml.  The  tumor  was  imbedded  in  the  intralipid  solution  of  0.7%  concentration,  and  the 
fitted  background  |ia  and  p,’ values  were  0.021  cm'1  and  6.83  cm'1,  respectively.  An  ultrasound 

B-scan  image  of  the  tumor  is  shown  in  Figure  4,  and  the  depth  of  the  tumor  evaluated  at  the 
center  of  the  mass  is  about  2.5  cm.  The  tumor  is  acoustically  heterogeneous  which  can  be  seen 
by  the  inhomogeneous  echo  pattern  in  the  image.  The  absorption  map  of  the  tumor  at  780  nm  is 
shown  in  Figure  5,  and  the  reconstructed  absorption  coefficient  is  0.17  cm'1.  To  evaluate  the 
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reconstructed  absorption  coefficients  at  different  target  volumes,  we  cut  the  tumor  by  half  step 
by  step  and  imaged  the  remaining  portions  of  the  tumor.  The  weights  were  2.36,  1.75,  1.12,  0.59 
grams,  respectively.  The  reconstructed  absorption  coefficients  were  0.18,  0.26,  0.27,  0.23  cm'1, 
respectively.  The  changes  in  |ia  are  primarily  related  to  the  heterogeneous  structure  of  the  tumor. 


Figure  4  .  Ultrasound  B-scan  of  the  excised  tumor  pointed  by  arrows.  Tumor  was  imbedded  in 
the  intralipid  of  approximately  0.7%  concentration  and  the  fitted  background  |iaand  (is  values 
are  0.021  cm'1  and  6.83  cm'1,  respectively. 

Key  Research  Accomplishments: 

•  refined  our  optical  imaging  system  and  implemented  imaging  software. 

•  optimized  our  combined  probe  by  studying  the  optical  and  ultrasound  sensor 
distributions. 

•  The  combined  method  has  been  evaluated  with  phantoms  and  excised  tumors  and  the 
improvement  of  the  method  compared  with  optical  or  ultrasound  alone  is  significant. 
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cm 

Figure  6.  Absorption  map  obtained  at  780  nm.  The  reconstructed  absorption  coefficient  is  0.17 
cm-1. 

Conclusions: 

We  have  refined  our  optical  imaging  system  and  implemented  imaging  software.  We  have  also 
optimized  our  combined  probe  by  studying  the  optical  and  ultrasound  sensor  distributions.  The 
combined  method  has  been  evaluated  with  phantoms  and  excised  tumors  and  the  improvement  of 
the  method  compared  with  optical  or  ultrasound  alone  is  significant. 

We  will  start  our  first  group  of  clinical  study  soon  after  we  recruit  several  patients  and  we  look 
forward  to  reporting  our  clinical  results  in  the  near  future. 
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Design  of  near-infrared  imaging  probe  with  the 
assistance  of  ultrasound  localization 


Quing  Zhu,  Nan  Guang  Chen,  Daqing  Piao,  Puyun  Guo,  and  XiaoHui  Ding 


A  total  of  364  optical  source- detector  pairs  were  deployed  uniformly  over  a  9  cm  X  9  cm  probe  area 
initially,  and  then  the  total  pairs  were  reduced  gradually  to  60  in  experimental  and  simulation  studies. 
For  each  source- detector  configuration,  three-dimensional  (3-D)  images  of  a  1-cm-diameter  absorber  of 
different  contrasts  were  reconstructed  from  the  measurements  made  with  a  frequency-domain  system. 
The  results  have  shown  that  more  than  160  source- detector  pairs  are  needed  to  reconstruct  the  absorp¬ 
tion  coefficient  to  within  60%  of  the  true  value  and  appropriate  spatial  and  contrast  resolution.  How¬ 
ever,  the  error  in  target  depth  estimated  from  3-D  images  was  more  than  1  cm  in  all  source- detector 
configurations.  With  the  a  priori  target  depth  information  provided  by  ultrasound,  the  accuracy  of  the 
reconstructed  absorption  coefficient  was  improved  by  15%  and  30%  on  average,  and  the  beam  width  was 
improved  by  24%  and  41%  on  average  for  high-  and  low-contrast  cases,  respectively.  The  speed  of 
reconstruction  was  improved  by  ten  times  on  average.  ©  2001  Optical  Society  of  America 
OCIS  codes:  170.0170,  170.3010,  170.5270,  170.7170,  170.3830.  ' 


1.  Introduction 

Recently,  optical  imaging  techniques  based  on  diffu¬ 
sive  near-infrared  (NIR)  light  have  been  employed  to 
obtain  interior  optical  properties  of  human  tissues.1"8 
Functional  imaging  with  NIR  light  has  the  potential 
to  detect  and  diagnose  diseases  or  cancers  through 
the  determination  of  hemoglobin  concentration,  blood 
02  saturation,  tissue  light  scattering,  water  concen¬ 
tration,  and  the  concentration  and  lifetime  of  exoge¬ 
nous  contrast  agents.  Optical  imaging  requires  that 
an  array  of  sources  and  detectors  be  distributed  di¬ 
rectly  or  coupled  through  optical  fibers  on  a  boundary 
surface.  Measurements  made  at  all  source- detector 
positions  can  be  used  in  tomographic  image  recon¬ 
struction  schemes  to  determine  optical  properties  of 
the  medium.  The  frequently  used  geometric  config¬ 
urations  of  sources  and  detectors  are  ring  arrays4’9-11 
and  planar  arrays.3*12-15  A  ring  array  consists  of 
multiple  sources  and  detectors  that  can  be  distrib¬ 
uted  uniformly  on  a  ring.  Optical  properties  of  the 
thin  tissue  slice  (two-dimensional  slice)  enclosed  by 
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the  ring  can  be  determined  from  all  measurements. 
A  planar  array  can  be  configured  with  either  trans¬ 
mission  or  reflection  geometries.  In  transmission 
geometry,  multiple  detectors  can  be  deployed  on  a 
planar  array,  and  multiple  sources  or  a  single  source 
can  be  deployed  on  an  opposite  plane  parallel  to  the 
detector  plane.  Optical  properties  of  the  three- 
dimensional  (3-D)  tissue  volume  between  the  source 
and  the  detector  planes  can  be  determined  from  all 
measurements.  In  reflection  geometry,  multiple 
sources  and  detectors  can  be  distributed  on  a  planar 
probe  that  can  be  hand-held.3*15  Optical  properties 
of  the  3-D  tissue  volume  at  slice  depths  below  the 
probe  can  be  determined  from  all  measurements. 
The  reflection  probe  configuration  is  desirable  for  the 
imaging  of  brain  and  breast  tissues. 

Although  many  researchers  in  the  field  have  con¬ 
structed  imaging  probes  using  reflection  geome¬ 
try,2*3*15  to  our  knowledge  the  required  total  number 
of  source- detector  pairs  over  a  given  probe  area 
needed  to  accurately  reconstruct  optical  properties 
and  localized  spatial  and  depth  distributions  has  not 
been  addressed  before.  In  this  paper  we  study  the 
relationship  between  the  total  number  of  source- 
detector  pairs  and  the  reconstructed  imaging  quality 
through  experimental  measurements.  Computer 
simulations  are  performed  to  assist  in  understanding 
the  experimental  results. 

Because  the  target  localization  from  diffusive 
waves  is  difficult,  our  group  and  others  have  intro¬ 
duced  use  of  a  priori  target  location  information  pro- 
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vided  by  ultrasound  to  improve  optical  imaging.1® 

In  this  paper  we  demonstrate  experimentally  that 
the  accurate  target  depth  information  can  signifi¬ 
cantly  improve  the  accuracy  of  the  reconstructed  ab¬ 
sorption  coefficient  and  the  reconstruction  speed  for 
any  optical  array  configuration. 

The  required  total  number  of  source- detector  pairs 
is  also  related  to  the  image  reconstruction  algorithms 
used.  In  this  paper  we  obtained  experimental  mea¬ 
surements  using  a  frequency-domain  system  with  the 
source  amplitude  modulated  at  140  MHz.  In  simu¬ 
lations,  forward  measurements  were  generated  by 
use  of  the  analytic  solution  of  a  photon  density  wave 
scattered  by  a  spherical  inhomogeneity  embedded  in 
a  semi-infinite  scattering  medium.19  In  both  exper¬ 
iments  and  simulations,  linear  perturbation  theory 
within  the  Born  approximation  was  used  to  relate 
optical  signals  at  the  probe  surface  to  absorption  vari¬ 
ations  in  each  volume  element  within  the  sample. 
The  total  least-squares  (TLS)  method20”22  was  used 
to  formulate  the  inverse  problem.  The  conjugate 
gradient  technique  was  employed  to  iteratively  solve 
the  inverse  problem.  Therefore  the  results  we  ob¬ 
tained  are  directly  relevant  to  the  probe  design  with 
reconstruction  algorithms  based  on  the  linear  pertur¬ 
bation  theory  and  can  bemsed  as  a  first-order  approx¬ 
imation  if  high-order  perturbations  are  employed  in 

image  reconstructions. 

This  paper  is  organized  as  follows.  In  Section  2 
we  describe  an  analytic: solution  used  to  generate 
simulated  forward  data;  the  Bom  approximation, 
and  the  TLS  method  for  image  reconstruction.  In 
Section  3  we  discuss  the  probe  geometry,  a  frequency- 
domain  NIR  system  used  to  acquire  the  experimental 
measurements,  and  an  ultrasound  subsystem  used  to 
acquire  the  target  depth  information,  computation 
procedufes  used  to  obtain  both  simulated  and  exper¬ 
imental  absorption  images.  In  Sections  4  and  5  we 
report  experimental  results  obtained  from  the  dense 
and  sparse .  arrays  with  and  without  a  priori  target 
depth  .information.  A  high:contrast  example  is 
given  in  Section  4,  and  a  low-contrast  case  is  given  in 
Section  5,  .  Simulations  are  performed  to  assist  in 

understanding  the' noise  on  the  image  reconstruction. 
Imaging  parameters,  evaluated  are  a  -6-dB  width  of 
the;  image  lobe,  the  reconstructed  maximum  value  of 
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the  image  artifact. .level.  ,  In.  Sections  6  and  7  we 
provide  a  discussion  and  summary,  respectively. 

*  '  i.  1  ~  -  j 

2.  Basic  PrinciplesJjt^jL  ; 

A.  Forward 

In  our  experimenfepforward  measurements  were 
made  with  a  frequency-domain  system  operating  at  a 
140-MHz  modulation -frequency.  In  computer  sim¬ 
ulations,  forward, -measurements  were  generated 
from  an  analytic  'solution  of  a photon  density  wave 
scattered  by  a  spherical  inhomogeneity . i9  When  the 
center  of  the  sphere';  coincides  with  the  origin  of 
spherical  coordinates  j.  the  solution  for  the  scattered 


iiiSllt 


Target 


Fig.  1.  Target,  source,  and  detector  configurations  for  a  semi¬ 
infinite  medium. 


photon  density  wave  Usc  outside  the  sphere  at  a  de¬ 
tector  position  r  =  (r,0,4>)  is  of  the  form 

UJjr,u)  =  2  {A^Mk^r)  +jnl(k°“tr)]YlJ*M, 

tm  (1) 


where  jt  and  are  spherical  Bessel  and  Neuman 
functions,  respectively;  Y/jm(0,4>)  are  the  spherical 
harmonics,  and  kout  =  [(“vp,aout  +  jco)/Dou  ]  is  the 
complex  wave  number  outside  the  sphere,  w  is  the 
angular  modulation  frequency  of  the  light  source, 


jxa°ut  is  the  absorption  coefficient  outside  the  sphere, 


and  Dout  is  the  photon  diffusion  coefficient  outside  the 
sphere  given  by  Dout  =  l/(3p,s'),  where  p.s'  is  the 
reduced  scattering  coefficient  outside  the  sphere. 
The  coefficients  A[>m,  determined  by  the  boundary 
conditions,  are 


Aijn  =  -UvSkoat/Dmt)h^\kmtrs)Y,,m*Ms) 
D°ulxji(y)ji'  (*)  -  Dinyji' ( y)ji(x ) 

X  D^xhPXxyfo)  -  D^hPKxWiy)  ’ 

where  x  =  kouta,  y  =  kina,  rs  =  (rs,0s,<t>s)  is  the  source 
position,  are  the  Hankel  functions  of  the  first 
kind,  and  j/  and  h\1]'  are  the  first  derivatives  of/,  and 
.  The  analytic  solution  has  the  important  adv an- 
tage  in  that  it  is  exact  to  all  orders  of  perturbation 
theory  and  thus  can  represent  accurate  measure¬ 
ments. 

We  generalized  the  above  analytic  solution  to  a 
semi-infinite  geometry  by  using  a  method  of  images 
with  extrapolated  boundary  conditions  (see  Fig.  1). 
A  type  I  boundary  condition  (zero  light  energy  density 
at  the  extrapolated  boundary)  is  used  to  derive  the 
scattered  wave  I/sc'.  To  calculate  the  Usc'  in  semi¬ 
infinite  geometry,  we  use  rsph+  =  (ro+,0o+,4>o+)  an^ 
rs  h_  =  (r0_,6o_,4>o-)  to  represent  the  centers  of  the 
sphere  and  the  image  sphere,  respectively.  The  vec¬ 
tors  r+'  =  r  -  rsph+  =  (r+',0+',4>+')  and  r_  =  r  - 
r  h_  =  (r_',0_',4)_')  are  therefore  pointing  to  the 
detector  position  r  =  (r,0,4>)  from  the  sphere  and  the 
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image  sphere,  respectively.  The  semi-infinite  solu¬ 
tion  of  Usc'  can  be  approximated  as 

U"'M  =  2  {Ai,m+[jt(k°^r+') 

l,m 

+  jnl(k°atr+')]YlJ*+',4>+')} 

~  2  {Ai,m-[ji(koulr+') 

l,m 

+  jni(koutr <!>+')} 

l,m 

+  jnl(k°atrY)]YljOY,  <>_')} 

-  2  {A^iMk^rJ) 

l,m 

+  jnl(koatrJ)]Yl^-',< (>_')},  (3) 

where 

At,m+  =  -0'uS^Vi?^1)(*ouV.+)  V(9.+,<l>.+) 
f  D°'itxjl(y)jl'(x)-Dinyjl’(y)ji(x) 


D^xhP'WMy)  -  ’ 

A;,m-  =  V(e.",4»._) 

Z?0Ut^(y)jV(x)  -  Dinyji  {y)Mx) 
Doatxh\1Y(x)My)  -  D^yh^ix)] i  (y) ' 

rs+  =  (rg+,es+,(|)s+)  and  rs“  =  (rs",es~,c|>s~)  are  the 
positions  of  the  source  and  the  image  source,  respec¬ 
tively. 

The  incident  photon  density  wave  at  the  detector 
position  r  has  the  following  form3: 


t/inc(r,co) 


S  [  exp(j^out|r  -  rs+|) 


exp(^out|r  -  r/|) 


The  total  photon  density  at  detector  r  is  a  super¬ 
position  of  its  incident  (homogeneous)  and  scattered 
(heterogeneous)  waves: 

U(  r,(o)  =  t/inc(r,co)  +  C/sc'(r,co).  (7) 

B.  Born  Approximation  for  Reconstruction 

The  Born  approximation  was  used  to  relate  Usc'  (r,co) 
measured  at  the  probe  surface  to  absorption  varia¬ 
tions  in  each  volume  element  within  the  sample.  In 
the  Born  approximation,  the  scattered  wave  that 
originated  from  a  source  at  rs  and  measured  at  de¬ 
tector  rd  can  be  related  to  the  medium  heterogeneity 
Afxa  (rv)  by 


U3  c'(rd,rs,a>)  =  G(rv,rrf,co)J7inc(rv,rs,a)) 


X  [vAfxa(rv)/Z)]dr v3, 


where  G{ rv,rrf,co)  is  the  Green  function  and  Ajma(rv)  = 
jjua(rv)  -  pa  is  the  medium  absorption  variation.11 
jla  is  the  average  value  of  the  medium  absorption 
coefficient.  By  breaking  the  medium  into  discrete 
voxels,  we  obtain  the  following  linear  equations: 

N 

Usc  (i*di,rsi, to)  —  G(rvj,rdi,w)  Ginc(iVj,rSi,o)) 

X  [v^a(rvj)/D]Arv\  (9) 

When  Wy  =  G(rvy,rdi,a))Uinc  (rty-,rsi,co)vA rv3/Z),  we  ob- 
tain  the  matrix  equation  of  Eq.  (9): 

[W]MJ^v{AfXa}iVXl  =  [UsdllMXl-  (10) 

The  realistic  constrains  on  Apa  are  (-a  X  back¬ 
ground  (LLa)  -  <  A[xa  <  1,  where  0  <  a  <1. 

The  above  constrains  ensure  that  the  reconstructed 
absorption  coefficient  jla  =  background  jjlg  +  Ap?  is 
positive  and  not  unrealistically  higher  than  unity. 
With  M  measurements  obtained  from  all  possible 
source- detector  pairs  in  the  planar  array,  we  can 
solve  N  unknowns  of  Ap,a  by  inverting  the  matrix  Eq. 
(10).  In  general,  the  perturbation  in  Eq.  (10)  is  un¬ 
derdetermined  (M  <  N)  and  ill-posed. 

When  the  target  depth  is  available  from  ultra¬ 
sound,  we  can  set  Ajjl^  of  a  nontarget  depth  equal  to 
zero.  This  implies  that  all  the  measured  perturba¬ 
tions  were  originated  from  the  particular  depth  that 
contained  the  target.  Because  the  number  of  un¬ 
knowns  was  reduced  significantly,  the  reconstruction 
converged  fast. 

C.  Total  Least-Squares  Solution 

To  solve  the  unknown  optical  properties  of  Eq.  (10), 
several  iterative  algorithms  have  been  used  in  the 
literature  including  the  regularized  least-squares 
method10  and  the  TLS  method.20*21  The  TLS  per¬ 
forms  better  than  the  regularized  least-squares 
method  when  the  measurement  data  are  subject  to 
noise  and  the  linear  operator  W  contains  errors. 
The  operator  errors  can  result  from  both  the  approx¬ 
imations  used  to  derive  the  linear  model  and  the 
numerical  errors  in  the  computation  of  the  operator. 
We  found  that  the  TLS  method  provides  more  accu¬ 
rate  reconstructed  optical  properties  than  the  regu¬ 
larized  least-squares  method,  so  we  adapted  the  TLS 
method  to  solve  the  inverse  problems.  It  has  been 
shown  by  Golub22  that  the  TLS  minimization  is 
equivalent  to  the  following  minimization  problem: 

.  IK/sd-wx||2 


where  X  represents  unknown  optical  properties. 
The  conjugate  gradient  technique  was  employed  to 
iteratively  solve  Eq.  (11). 

3.  Methods 

A.  Probe  Design  and  Imaging  Geometry 

There  are  two  basic  requirements  to  guide  the  design 
of  the  NIR  probe.  First,  all  source- detector  separa- 


3290 


APPLIED  OPTICS  /  Vol.  40,  No.  19  /  1  July  2001 


Fig  2  Configuration  of  a  dense  array  with  28  optical  sources  and 
13  detectors  as  well  as  six  ultrasound  transducers.  Large  black 
circles  are  optical  detectors,  gray  circles  are  optical  sources,  and 
small  white  circles  are  ultrasound  transducers.  A  1-cm-diameter 
spherical  target  was  located  at  various  depths  in  simulations  and 
experiments. 


Fig.  3.  Schematic  of  a  single-channel  optical  data-acquisition  sys¬ 
tem  A 140  02-MHz  oscillator  is  used  to  drive  the  laser  diode  (78U 
nm)'  that  delivers  the  light  to  the  medium  through  the  fiber.  The 
detected  signals  are  amplified  and  mixed  with  signals  from  a  140- 
MHz  oscillator.  The  heterodyned  20-kHz  signals  are  amplified, 
filtered,  and  digitized.  The  signals  from  two  oscillators  are  also 
directly  mixed  to  provide  reference  signals.  The  amplitude  and 
phase  of  the  waveform  received  through  the  medium  are  calculated 
from  signals  measured  through  the  medium  and  the  reference. 
PMT  nhntomultiolier  tube. 


tions  have  to  be  as  large  as  1  cm  so  that  diffusion 
theory  is  a  valid  approximation  for  image  reconstruc¬ 
tion.  Second,  because  the  depth  of  a  photon  path  is 
measured  approximately  one  third  to  one  hall  the 
source- detector  separation,  the  distribution  ol 
source-detector  distances  should  be  from  approxi¬ 
mately  1  to  10  cm  to  effectively  probe  the  depth  from 
approximately  0.5  to  4  cm.  On  the  basis  of  these 
requirements,  we  deployed  a  total  of  28  sources  and 
13  detectors  over  a  probe  area  of  9  cm  x  9  cm  (see  1  ig. 
2)  The  minimum  source- detector  separation  in  the 

configuration  is  1.4  cm  and  the  maximum  is  10.0  cm. 
We  call  this  array  a  filled  or  dense  array  (a  term 
adapted  from  ultrasound  array  design).  The  9  cm  X 
9  cm  X  4  cm  imaging  volume  was  discretized  into 
voxels  of  size  0.5  cm  X  0.5  cm  X  1  cm;  therefore  a  total 
of  four  layers  in  depth  was  obtained.  The  target  was 
a  1-cm-diameter  sphere  located  at  different  locations. 
Because  one  of  the  objectives  of  this  study  was  to 
evaluate  the  target  depth  distribution,  the  centers  of 
the  four  layers  in  depth  were  adapted  to  the  target 
depth.  For  example,  if  the  target  depth  was  z  -  3 
cm,  the  centers  of  the  four  layers  were  chosen  as  1,  2, 
3,  and  4  cm,  respectively. 

The  ultrasound  transducers  shown  in  Fig.  I  were 
deployed  simultaneously  on  the  same  probe,  lhe 
diameter  of  each  ultrasound  transducer  is  1.5  mm 
and  the  spacing  between  the  transducers  is  4  mm, 
except  the  two  located  closer  to  the  optical  detector  in 
the  middle.  The  spacing  between  these  two  trans¬ 
ducers  is  8  mm.  Because  this  study  requires  accu¬ 
rate  target  location  as  a  reference  to  compare  with 
the  reconstructed  absorption  image  location,  six 
transducers  are  used  to  guide  the  spatial  positioning 
of  a  target.  The  target  is  centered  when  the  two 


middle  ultrasound  transducers  receive  the  strongest 
signals.  The  target  depth  is  determined  from  re¬ 
turned  pulse-echo  signals.  In  this  study  we  do  not 
intend  to  provide  ultrasound  images  of  the  target 
with  such  a  sparse  ultrasound  array,  but  we  demon¬ 
strate  the  feasibility  of  using  a  priori  depth  informa¬ 
tion  to  improve  optical  reconstruction. 

In  regard  to  the  image  voxel  size,  there  is  a  trade¬ 
off  between  the  accurate  estimation  of  the  weight 
matrix  W  and  the  voxel  size.  Because  Wy-  is  a  dis¬ 
crete  approximation  of  the  integral 

’G(rv,rd,ai)f/inc(r„,rs,a))(v/D)drv3, 

Vu 


it  is  more  accurate  when  the  voxel  size  is  smaller. 
However,  the  total  number  of  reconstructed  un¬ 
knowns  will  increase  dramatically  with  the  decreas¬ 
ing  voxel  size.  Because  the  rank  of  the  matrix  W  is 
less  than  or  equal  to  the  total  number  of  measure¬ 
ments  [Eq.  (10)],  a  further  decrease  m  voxel  size  will 
not  add  more  independent  information  to  the  weight 
matrix.  We  found  that  a  0.5  cm  X  0.5  cm  X  1  cm 
voxel  size  is  a  good  compromise.  Therefore  we  used 
this  voxel  size  in  image  reconstructions  reported  in 
this  paper. 

B.  Experimental  System 

We  constructed  a  NIR  frequency-domain  system  and 
the  block  diagram  of  the  system  is  shown  in  Fig.  3.  Un 
the  source  side,  a  140.000-MHz  sine-wave  oscillator 
was  used  to  modulate  the  output  of  a  780-nm  diode 
laser  that  was  housed  in  an  optical  coupler  (OZ,  (Jp- 
tics  Inc.).  The  output  of  the  diode  was  coupled  to  the 
turbid  medium  through  a  single  200- pm  multimode 
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optic  fiber.  On  the  reception  side,  an  optical  fiber  of 
3  mm  in  diameter  was  used  to  couple  the  detected 
light  to  a  photomultiplier  tube  detector.  The  output 
of  the  photomultiplier  tube  was  amplified  and  then 
mixed  with  a  local  oscillator  at  a  frequency  of  140.020 
MHz.  The  heterodyned  signal  at  20  kHz  after  the 
mixer  was  further  amplified  and  filtered  by  a  band¬ 
pass  filter.  The  outputs  of  two  oscillators  (140.000- 
and  140.020-MHz  signals)  were  directly  mixed  to 
produce  20-kHz  reference  signals.  Both  signals 
were  sampled  simultaneously  by  a  dual-channel  1.25- 
MHz  analog-to-digital  converter  (A/D)  board.  The 
Hilbert  transform  was  performed  on  both  sampled 
and  reference  waveforms.  The  amplitude  of  the  Hil¬ 
bert  transform  of  the  sampled  waveform  corresponds 
to  the  measured  amplitude,  and  the  phase  difference 
between  the  phases  of  the  Hilbert  transforms  of  the 
sampled  and  reference  waveforms  corresponds  to  the 
measured  phase. 

A  black  probe  with  holes  shown  in  Fig.  2  was  used 
to  emulate  the  semi-infinite  boundary  condition. 
Two  3-D  positioners  were  moved  independently  to 
position  the  source  and  detector  fibers  at  the  desired 
spatial  locations  within  the  9  cm  X  9  cm  area. 

A  challenge  in  the  reflection  NIR  probe  design  is  to 
preserve  a  huge  dynamic  range  in  received  signals. 
For  example,  the  amplitude  at  a  1-cm  source- 
detector  separation  measured  from  0.6%  Intralipid  in 
reflection  mode  is  approximately  84  dB  larger  than 
that  at  a  9-cm  separation.  So  the  signals  can  be 
saturated  when  they  are  measured  from  closer 
source- detector  pairs,  but  they  may  be  too  low  at 
more  distant  source-detector  pairs.  The  problem 
can  be  overcome  by  means  of  controlling  the  light 
illumination.  At  least  two  illumination  conditions 
need  to  be  used:  a  low  source  level  for  closer  source- 
detector  pairs  and  a  high  level  for  distance  source- 
detector  pairs.  In  our  system,  a  30-dB  attenuator 
connected  to  the  140-MHz  oscillator  was  switched  on 
and  off  to  provide  two  different  source  levels  and  thus 
to  preserve  the  dynamic  range.  Figure  4(a)  shows  a 
plot  of  the  measured  log  [p2t/sd(p)]  versus  the  source- 
detector  separation  p,  and  Fig.  4(b)  shows  the  plot  of 
the  measured  phase  versus  the  source- detector 
separation.  The  Intralipid  concentration  was 
0.6%,  which  corresponded  to  \±s'  =  6  cm-1.  Theo¬ 
retically  both  log  [p2C/sd(p)]  and  phase  are  linearly 
related  to  the  source- detector  separation  because 
of  the  semi-infinite  boundary  condition  used,3  and 
experimental  measurements  shown  in  Fig.  4  vali¬ 
date  that  they  are  linearly  related  to  the  source- 
detector  separation. 

The  ultrasound  system  consists  of  six  transducers 
(see  Fig.  5),  a  pulser  (Panametrics  Inc.),  an  A/D  con¬ 
verter,  and  a  multiplexer.  The  pulser  provided  a 
high-voltage  pulse  of  6-MHz  central  frequency  to  drive 
each  selected  transducer.  The  returned  signals  were 
received  by  the  same  transducer,  amplified  by  the  re¬ 
ceiving  circuit  inside  the  pulser,  and  sampled  by  the 
A/D  converter  with  a  100-MHz  sampling  frequency. 


Source- Detector  Separation  (cm) 


Source-Detector  Separation  (cm) 

Fig.  4.  Calibration  curves,  (a)  log  [p2U(p)]  versus  source- 
detector  separation,  (b)  Phase  versus  source- detector  separa¬ 
tion. 


C.  Computation  Procedures 

L  Computation  Procedures  of  Experimental  Data 
To  study  the  relationship  between  the  total  number 
of  source- detector  pairs  and  distributions  of  recon¬ 
structed  optical  absorption  coefficients,  we  started 
from  the  dense  array  with  a  total  of  28  sources  and  13 
detectors  (see  Fig.  2)  and  gradually  reduced  this 
number  to  generate  sparse  arrays  with  24  X  13  (24 
sources  and  13  detectors),  20  X  13,  28  X  9,  24  X  9, 
16  X  13,  20  X  9, 12  X  13, 16  X  9,  28  X  5,  24  X  5,  12  X 
9,  20  X  5,  16  X  5,  and  12  X  5  source- detector  pairs, 
respectively.  Each  sparse  array  was  a  subset  of  the 


Ultrasound 


Fig.  5.  Ultrasound  data-acquisition  system.  The  pulser  is  used 
to  generate  high-voltage  pulses  that  are  used  to  excite  the  selected 
ultrasound  transducer.  The  returned  signals  are  received  by  the 
selected  transducer  and  are  sampled  by  the  A/D  converter. 
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dense  array,  and  its  probe  area  was  the  same  as  the 
dense  array.  For  each  sparse  array  configuration, 
we  compared  the  reconstructed  optical  imaging  pa¬ 
rameters  measured  from  the  dense  array  with  those 
from  the  sparse  arrays.  The  parameters  include  the 
maximum  values  of  reconstructed  |ia  at  different  lay¬ 
ers  and  their  spatial  locations,  spatial  resolution  and 
artifact  level  of  the  |i0  image,  and  target  depth  dis¬ 
tribution.  Targets  of  different  absorption  contrasts 
were  located  at  different  positions  inside  the  In¬ 
tralipid.  For  each  target  case,  one  set  of  measure¬ 
ments  with  the  dense  array  was  obtained,  and 
subsets  of  the  measurements  were  used  as  sparse 
array  measurements.  In  all  experiments,  the  back- 
ground  Intralipid  concentration  was  approximately 
0.6%,  and  |xa.'  was  experimentally  determined  from 
curve  fitting  results.  Currently,  we  did  not  recon¬ 
struct  target  p.s',  and  we  used  the  common  \±J  for 
both  the  background  and  the  target. 

The  total  number  of  iterations  or  stopping  criterion 
was  difficult  to  determine  for  experimental  data. 
Ideally,  the  iteration  should  stop  when  the  object 
function  [see  Eq.  (11)]  or  the  error  performance  sur¬ 
face  reaches  the  noise  floor.  However,  the  system 
noise,  particularly  coherent  noise,  was  difficult  to  es¬ 
timate.  In  general,  we  found  that  the  reconstructed 
values  were  closer  to  true  values  when  the  object 
function  reached  approximately  5-15%  of  the  initial 
value  (total  energy  in  the  measurements).  How¬ 
ever,  this  criterion  was  applicable  only  to  reconstruc¬ 
tions  with  total  source- detector  pairs  closer  to  the 
dense  array  case.  Therefore  we  used  this  criterion 
for  the  dense  array  reconstruction  and  used  the  same 
iteration  number  obtained  from  the  dense  array  for 
the  sparse  arrays.  Thus  the  iteration  number  is 
normalized  to  the  dense  array  case. 


2.  Computation  Procedures  of  Simulation 

Simulations  were  performed  to  assist  the  under¬ 
standing  of  the  random  noise  on  the  reconstructed 
absorption  coefficient.  In  simulations,  Gaussian 
noise  with  different  standard  deviations  proportional 
to  the  average  value  of  each  forward  data  set  was 
added  to  the  forward  measurements.  Typically, 
0.5%,  1.5%,  and  2.0%  of  the  average  value  of  each 
forward  data  set  were  used  as  standard  deviations  to 
generate  noise.  In  simulations,  the  target  |x0  was 
changed  to  different  contrast  values,  and  target  p,s 
was  kept  the  same  as  the  background.  The  back¬ 
ground  |i.n  and  |V  were  0.02  and  6  cm  \  respectively. 

The  stopping  criterion  used  in  the  simulations  was 
based  on  the  noise  level  of  the  object  function.  Con¬ 
sidering  that  the  object  function  fluctuates  within  one 
standard  deviation  ct  around  the  mean  E  when  the 
iteration  number  is  large,  we  can  use  E  +  cr  as  a 
stopping  criterion,  i.e.,  the  iteration  will  stop  it  the 
object  function  is  less  than  E  +  cr.  When  the  linear 
perturbation  is  assumed,  E  can  be  approximated  as 
Sfn(j)2 *  and  cr  as  2f[2 n(y)4]1/2,  where  N  is  the  total 
number  of  source- detector  pairs  and  n(j)  is  the  gen¬ 
erated  random  noise  with  a  standard  deviation  pro 


portional  to  the  specified  percentage  of  the  mean  of 
the  forward  data  set. 


D.  Testing  Targets 

Spherical  testing  targets  of  ~1  cm  in  diameter  were 
made  of  acrylamide  gel.16  The  acrylamide  powder 
was  dissolved  in  distilled  water,  and  20%  concentra¬ 
tion  of  Intralipid  was  added  to  the  acrylamide  solu¬ 
tion  to  dilute  the  solution  to  a  0.6%  Intralipid 
concentration  (|xs'  =  6  cm  1).  India  ink  was  added 
to  the  solution  to  produce  target  |xa  of  different  val¬ 
ues.  Acoustic  scattering  particles  of  200  p-m  in  di¬ 
ameter  were  added  to  the  solution  before 
polymerization.  Components  of  ammonium  persul¬ 
fate  and  tetramethylethylenediamine  (known  as  TE- 
MED)  were  added  to  the  solution  to  produce 
polymerization. 


4.  Results  of  a  High-Contrast  Target  Case 

A.  Experimental  Results  of  a  Dense  Array 
Figure  6(a)  is  an  experimental  image  of  a  high- 
contrast  target  (pa  =  0.25  cm  l)  located  in  the  In¬ 
tralipid  background  (ps'  =  6  cm  ).  The  target  was 
a  1-cm-diameter  sphere  and  its  center  was  located  at 
(x  =  0,  y  =  0,  z  =  3.0  cm),  where  x  and  y  were  the 
spatial  coordinates  and  z  was  the  propagation  depth. 
The  target  depth  was  well  controlled  by  use  of  ultra¬ 
sound  pulse-echo  signals,  and  the  error  was  less  than 
1  mm.  The  3-D  images  were  reconstructed  from  the 
measurements  made  with  the  dense  array,  and  the 
image  shown  was  obtained  at  target  layer  3.  The 
centers  of  the  imaging  voxels  in  z  are  1,2,3,  and  4  cm 
for  layers  1,  2,  3,  and  4,  respectively.  The  measured 
maximum  value  of  the  image  lobe  [fla(max)] 0.233 
which  was  a  close  estimate  of  the  target 
The  measured  spatial  location  of  |ia(max)  was  {x  =  0.5 
cm,  y  =  0.0  cm),  which  agreed  reasonably  well  with 
the  true  target  location.  The  spatial  resolution  can 
be  estimated  from  the  -6-dB  contour  plot  of  Fig.  6(a), 
which  is  shown  in  Fig.  6(b) .  The  outer  contour  is  -6 
dB  from  the  |i0(max)>  and  the  contour  spacing  is  1  d  . 
The  width  of  the  image  lobe  measured  at  the  -b 
dB-level  corresponds  to  a  full  width  at  half-maximum 
(FWHM),  which  is  commonly  used  to  estimate  reso¬ 
lution.  The  measured  widths  of  longer  and  shorter 
axes  were  1.01  and  1.60  cm,  respectively,  and  the 
geometric  mean  was  1.27  cm,  which  was  used  to  rep¬ 
resent  the  -6-dB  beam  width.  The  contrast  resolu¬ 
tion  can  be  estimated  from  the  peak  artifact  level, 
and  no  artifact  was  observed  in  the  image.  The  tar¬ 
get  depth  can  be  assessed  from  the  images  obtained 
from  other  nontarget  layers.  Figure  6(c)  is  the  im¬ 
age  obtained  at  nontarget  layer  4,  and  an  image  lobe 
of  lL  ,  —  0.138  cm'1  was  observed.  The  spatial 

location* of  |ia(max)  was  (x  =  0.0,  y  =  0.0),  which 
agreed  well  with  the  true  target  location.  Mo  dis¬ 
tinct  lobes  were  observed  at  nontarget  layers  1  and  2, 
which  indicates  that  the  error  in  the  target  depth 
estimated  from  3-D  images  was  approximately  1  cm. 
Because  the  error  in  the  true  target  depth  was  less 
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Fig.  6.  Experimental  3-D  images  of  fxa  reconstructed  with  a  total 
of  28  x  13  =  364  source- detector  pairs  at  2437  iterations.  The 
target  (fxa  =  0.25  cm-1)  was  located  at  {x  =  0,  y  ~  0,  z  -  3.0  cm) 
inside  the  Intralipid  background,  (a)  Reconstructed  jia  at  target 
layer  3.  The  horizontal  axes  represent  spatial  *  andy  coordinates 
in  centimeters,  and  the  vertical  axis  is  the  j la.  The  measured 
maximum  value  of  the  image  lobe  [fia(max)]  was  0.233  cm"1,  and  its 
location  was  (;e  —  0.5,  y  =  0.0).  No  image  artifacts  were  observed, 
(b)  -6-dB  contour  plot  of  (a).  The  outer  contour  is  -6  dB  from  the 
M'a(mnx))  and  the  contour  spacing  is  1  dB.  The  measured  -6-dB 
beam  width  was  1.27  cm.  (c)  Reconstructed  }ia  at  nontarget  layer 
4.  An  image  lobe  of  strength  0.138  cm-1  and  spatial  location  of 
( x  —  0.0,  y  =  0.0)  was  observed. 


than  1  mm,  this  1-cm  error  was  due  largely  to  the 
depth  uncertainty  of  diffusive  waves. 

B.  Simulation  Results  of  a  Dense  Array 

Our  simulations  support  the  experimental  results. 
Figure  7  shows  simulation  results  obtained  with  the 
dense  array.  A  simulated  1-cm-diameter  absorber 
(fxa  =  0.25  cm  *)  was  located  at  (x  =  0,  y  =  0,  z  =  3 
cm)  inside  a  homogeneous  scattering  background  (jx</ 
=  6  cm  1).  We  added  0.5%  Gaussian  noise  to  the 
forward  data  generated  from  the  analytic  solution. 
Images  obtained  at  nontarget  layer  2,  target  layer  3, 
and  nontarget  layer  4  are  shown  in  Fig.  7(a),  Fig.  (b), 
and  7(c),  respectively.  No  target  was  found  at  layer 


Fig.  7.  Simulated  3-D  images  of  |la  reconstructed  with  a  total  of 
28  X  13  =  364  source-detector  pairs.  The  target  (jul0  =  0.25  cm-1) 
was  located  at  (x  —  0,  y  =  0,  z  =  3.0  cm)  inside  the  Intralipid 
background,  (a)  Reconstructed  |ia  at  nontarget  layer  2  (simula¬ 
tion,  0.5%  noise).  No  image  lobe  was  observed,  (b)  Recon¬ 
structed  |ia  at  target  layer  3  (simulation,  0.5%  noise).  The  target 
of  strength  |!a(max)  -  0.248  cm-1  and  the  spatial  location  (0.0,  0.0) 
was  observed,  (c)  Reconstructed  jia  at  nontarget  layer  4  (simula¬ 
tion,  0.5%  noise).  The  target  of  strength  0.190  cm"1  and  location 
(x  =  0.0,  y  =  0.0)  was  observed,  (d).  Reconstructed  j la  at  non¬ 
target  layer  2  with  1.0%  noise.  The  target  of  0.028  cm'1  was 
observed.  Note  that  the  scale  of  (d)  is  different  from  (a)-(c). 


2,  and  the  target  of  strengths  0.248  and  0.190  cm-1 
appeared  at  layers  3  and  4,  respectively.  However, 
when  the  noise  level  in  the  forward  data  was  in¬ 
creased  to  1.0%,  the  target  of  strength  jia(max)  =  0.028 
cm  appeared  at  nontarget  layer  2  [Fig.  7(d)]  as  well 
as  target  layer  3  [|la(max)  —  0.163  cm-1]  and  nontarget 
layer  4  [|ia(max)  =  0.111  cm-1].  This  suggests  that 
measurement  noise  is  an  important  parameter  to  af¬ 
fect  the  target  depth  estimate. 

C.  Experimental  Results  of  a  Dense  Array  with 
Ultrasound  Assistance 

From  ultrasound  we  obtained  the  target  depth  as 
well  as  the  target  boundary  information.  Figure 
8(a)  shows  the  received  pulse-echo  signals  (from  a 
depth  of  2.3— 3.5  cm)  obtained  from  the  six  ultrasound 
transducers  (see  Fig.  2).  The  spatial  dimension  cov¬ 
ered  by  the  six  transducers  is  2.4  cm.  The  front 
surface  of  the  target  was  indicated  clearly  by  the 
returned  pulses  shown  by  two  arrows,  and  the  back 
surface  was  seen  through  the  reflection  of  a  soft  plas¬ 
tic  plate.  The  reflected  signals  of  the  back  surface 
are  also  shown  by  arrows.  The  plate  was  used  to 
hold  the  target  and  was  transparent  to  light.  The 
measured  depth  of  the  target  front  surface  was  2.44 
cm  and  the  back  surface  was  3.43  cm.  The  distance 
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Fig.  8.;  (a)  Ultrasound  j)ulse-echo  signals  or  A-scan  lines  obtained 
frorii  six  transducers.  The  abscissa  is  the  propagation  depth  in 
millimeters:  From  reflected  signals,  the  measured  depth  of  the 
target  front  surface  is  2.44  cm,  and  the  back  surface  is  3.43  cm. 
The  center  of  the  target  is  "3  cm.  The  total  length  of  the  signal 
corresponds  to  1.2  cm  in  depth,  and  the  measured  distance  be¬ 
tween  the  front  and  the  back  surfaces  is  0.993  cm.  The  spatial 

dimension  covered  by  the  transducers  is  2.4  cm.  (b)  An  image  of 
the  high-contrast  target  reconstructed  at  a  target  layer  when  we 
used  only  a  priori  target  depth  information  provided  by  ultra- 
sbund.'  The  reconstructed  (lacAax)  reached  0.245  cm  1  at  216  it- 
erations;’ ' '■ v 

between  the  peaks  of  front  reflection  and  backreflec- 
tiofi  was  0.993.  cm,  which  corresponded  to  the  target 
size;  Wi^thei^ssistance'  of  target  depth,  we  recon¬ 
structed’  thie  absorption  image  at  the  target  layer 
only.'  1  Figure’.8(b)  is  the  reconstructed  |ia  obtained 
from  the  dense; aiTay  measurements.  A  total  of  2 16 
iterations  ^refused  to  obtain  |ia(max)  =  0.245  cm-1, 
and  the  reconstnictioh  was  approximately  ten  times 
faster  the  depth  information.  The 

spatial "reSolution^as  approximately  the  same  as 
t'ig**  rbreato  width  was  1.31  cm. 

The  cbntrasf'i’esblution  was  5  dB  poorer  because  the 
measurement  lumped  to  single-layer  re- 

distributed  to  four  lay- 


. rr  rCr  l  •  V.  •  -)  i  . 

D.  Array 

The  imdgi^g’^d0^y:-  of  sparse  arrays  decreased.  Fig¬ 
ure  ^(aX^l^^xperiinental  image  at  target  layer  3 
obtained  X  5  sparse  array.  The  sparse 

array  measurements  used  were  a  subset  of  the  dense 


Fig.  9.  Experimental  images  at  target  layer  3  reconstructed  with 
a  total  of  16  X  5  —  80  source- detector  pairs,  (a)  Reconstructed  |la 
at  target  layer  3  with  2478  iterations.  The  measured  m^max) was 
0  107  cm"1,  which  was  43%  of  the  true  value,  and  the  spatial 
location  of  |ia(max)  was  (*  =  0.5,  y  »  -0.5).  The  measured  peak 
image  artifact  level  was  -10  dB  below  the  peak  of  the  mam  image 
iobe.  (b)  —  12-dB  contour  plot  of  (a).  The  outer  contour  is  -12 
dB,  and  the  contour  spacing  is  2  dB.  The  measured  —  6*dB  beam 
width  was  2.55  cm,  which  was  200%  broader  than  that  of  the  dense 
array,  (c)  Reconstructed  (la  at  target  layer  3  with  10,000  itera¬ 
tions.  The  measured  |l0(max)  reached  0.264  cm  \  and  the  peak 
artifact  level  was  increased  by  2  dB  as  well,  (d)  Reconstructed  }la 
at  the  target  layer  with  only  a  priori  target  depth  information 
provided  by  ultrasound.  |ia(max)  =  0.173  cm  1  at  216  iterations. 


array  measurements.  The  measured  |ia(niax)  was 
0.107  cm"1,  which  was  43%  of  the  true  value.  The 
-6-dB  beam  width  was  2.55  cm,  which  was  200% 
broader  than  that  of  the  dense  array.  Edge  artifacts 
were  observed  and  are  best  seen  from  Fig.  9(b),  which 
is  the  -12-dB  contour  plot  of  Fig.  9(a).  The  outer 
contour  is  — 12  dB,  and  the  spacing  is  2  dB.  The  peak 
artifact  level  is  -10  dB  from  the  |ia(max)*  The  recon¬ 
structed  jia  can  be  increased  if  the  iteration  is  signifi¬ 
cantly  increased.  The  iteration  number  used  to 
obtain  Fig.  9(a)  was  2478,  which  was  the  same  as  that 
used  to  obtain  Fig.  6.  When  the  iteration  number  was 
increased  to  10,000  for  the  sparse  array  case,  the  re¬ 
constructed  fia(max)  reached  0.264  cm  1,  which  was 
close  to  the  true  target  \xa.  However,  the  image  arti¬ 
fact  level  was  increased  too  [see  Fig.  9(c)],  and  the  ratio 
of  the  peak  artifact  to  |ia(max)  was  2  dB  higher  than 
that  shown  in  Fig.  9(a).  In  addition,  the  background 
noise  fluctuation  of  nontarget  layers  1  and  2  was  in¬ 
creased.  The  noise  fluctuation  can  be  estimated  from 
the  standard  deviations  of  reconstructed  jla  at  nontar- 
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get  la3^ers  1  and  2.  At  2478  iterations,  the  averages 
and  the  standard  deviations  of  jla  were  0  0233 
(±0.0028)  and  0.0202  cm'1  (±0.0022  cm’1)  for  non¬ 
target  layers  1  and  2,  respectively,  and  these  values 
were  0.0244  (±0.0078)  and  0.0208  cm’1  (±0.0073 
cm  x)  at  10,000  iterations. 

The  measured  maximum  values  of  the  image  lobes 
or  target  strengths  at  the  target  layer  and  nontarget 
layer  4  continuously  grew  with  each  iteration,  even 
though  the  reconstructed  values  at  the  target  layer 
were  close  to  the  true  value.  This  problem  was  men¬ 
tioned  in  the  literature,11  but  was  not  explained  well. 
It  is  largely  related  to  use  of  nonlinear  constrains  on 
Apa  in  Eq.  (10),  particularly  the  choice  of  a.  When  a 
is  close  to  1,  the  reconstruction  converges  fast,  and 
the  target  strength  increases  little  after  a  certain 
number  of  iterations.  When  a  is  close  to  0,  the  re¬ 
construction  converges  slowly,  and  the’  target 
strength  grows  continuously.  However,  when  the 
measurement  signal-to-noise  ratio  (SNR)  is  not  high, 
for  example,  in  sparse  array  or  low-contrast  cases’ 
the  choice  of  a  **  1  can  cause  unstable  reconstruction! 
m  some  cases,  the  reconstructed  images  can  jump 
from  one  set  of  (xa  to  another,  which  causes  the  object 
junction  to  increase  suddenly  and  reduce  again. 
This  is  related  to  the  underdetermined  nature  of  Eq. 
(10),  i.e.,  the  unknowns  are  far  more  than  the  mea¬ 
surements.  In  some  cases  the  reconstructed  images 
have  multiple  lobes  of  similar  strengths,  which  indi¬ 
cate  that  the  reconstruction  does  riot  converge  at  all. 
In  all  cases,  the  image  background  fluctuations  were 
large  compared  with  the  fluctuations  when  a  smaller 
a  was  used.  We  found  that  a  between  0.1  and  0.4 
can  provide  stable  reconstruction,  and  we  used  a  =» 
0.1  for  all  experiments. 

Another  factor  that  accounts  for  the  slow  increase 
in  the  reconstructed  value  is  use  of  linear  perturba- 
tion  to  approximate  the  measurements  that  contain 
all  higher-order  perturbations.  The  minimization 
procedure  [Eq.  (11)]  blindly  minimizes  the  difference 
between  the  measurements  and  their  linear  approx¬ 
imation  WAjxa  and  therefore  reconstructs  higher  and 
higher  A|x0  if  the  iteration  continues. 

The  target  depth  estimate  was  poorer  than  that  of 
the  dense  array  because  of  the  lower  SNR  of  the 
sparse  array  measurements.  Similar  to  the  dense 
array  case,  a  target  of  strength  |ia(max)  =  0.072  cm"1 
and  location  ( x  =  0.5,  y  =  0.0)  was  observed  at  non¬ 
target  layer  4.  In  addition,  a  target  of  strength 
(max)  ~  0.0348  cm  1  and  location  (x  =  0.5,  y  = 
-0.5)  was  observed  at  nontarget  layer  2.  However, 
the  target  mass,  which  was  approximately  the  vol¬ 
ume  underneath  the  image  lobe,  was  much  smaller 
than  that  obtained  at  layers  3  and  4,  and  the  target  at 
layer  2  was  buried  in  the  background  noise. 

Figure  9(d)  is  the  reconstructed  |ia  at  the  target 
layer  from  only  the  sparse  array  measurements.  A 
total  of  216  iteration  steps  were  used  to  obtain 
Fa(max)  =  0.173  cm  1,  and  the  reconstruction  was 
approximately  50  times  faster  than  that  without  the 
depth  information.  The  measured  —  6-dB  beam 
width  was  1.49  cm,  which  was  approximately  the 
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Fig.  10.  |l0(max)  versus  the  total  number  of  source— detector  pairs. 
The  center  of  the  target  (p.a  =  0.25  cm-1)  was  located  at  (x  =  0.0, 
y  =  0.0, 2  =  3.0  cm)  in  computer  simulations  and  experiments,  (a) 
Curves  were  obtained  at  the  target  layer.  Two  dashed  curves 
(upper  and  lower)  are  the  curve-fitting  results  of  simulation  data 
points  obtained  with  0.5%  and  2.0%  noise  added  to  the  forward 
data,  respectively.  The  experimental  data  are  plotted  with  cir¬ 
cles,  and  the  dashed  curve  in  the  middle  is  the  fitting  result  of  the 
experimental  points,  (b)  The  measured  target  strength  (circles) 
and  the  curve-fitting  result  (lower  curve).  The  measured  target 
strength  (stars)  was  reconstructed  at  the  target  layer  only  by  use 
of  a  priori  depth  information  and  the  curve  fitting  result  (upper 
curve). 


same  as  Fig.  9(c)  but  was  improved  42%  from  Fig. 
9(a).  The  contrast  resolution  was  2  dB  worse  for  the 
same  reason  discussed  above. 

E.  Simulation  and  Experimental  Results  of  Reconstructed 
fla(max)  versus  Total  Number  of  Source-Detector  Pairs 

To  understand  the  effects  of  random  noise  on  the 
performance  of  the  reconstruction,  we  performed 
simulations  for  each  array  configuration.  Gauss¬ 
ian  noise  of  0.5%,  1.0%,  1.5%,  and  2.0%  were  added 
to  each  forward  data  set  generated  from  the  ana¬ 
lytic  solution  [see  Eq.  (3)].  The  center  of  a  simu¬ 
lated  1-cm-diameter  spherical  target  (pa  =  0  25 


cm-1)  was  located  at  {x  =  0,  y  =  0,  z  =  3.0  cm). 
Reconstructed  images  at  different  noise  levels  were 
obtained  for  each  array  configuration,  and  the  peak 
values  of  image  lobes  [Aa(max)]  target  layer  3  were 
measured.  Figure  10(a)  shows  simulation  and  ex¬ 
perimental  results  of  reconstructed  M-at max)  versus 
the  total  number  of  source- detector  pairs.  Two 
dashed  curves  are  the  fitting  results  of  simulated 
data  points  with  0.5%  (upper)  and  2.0%  (lower) 
noise.  The  dashed  curve  in  the  middle  is  the  fitting 
result  of  experimental  points  plotted  with  circles. 
Second-order  polynomials  were  used  for  all  curve 
fittings.  The  reduction  of  the  reconstructed  p.a(max) 
was  significant  when  the  noise  level  went  up  in  the 
simulated  data.  Because  the  SNR  of  experimental 
data  was  decreased  when  the  total  number  of 
source— detector  pairs  was  reduced,  the  data  were 
scattered  around  the  0.5%  noise  curve  when  the 
total  pairs  were  large  and  were  distributed  around 
the  2.0%  noise  curve  (values  were  60%  less  than 
that  obtained  from  the  dense  array)  when  the  total 
pairs  were  reduced  to  less  than  140.  However,  be¬ 
cause  our  experimental  system  has  both  coherent 
and  random  noise,  simulations  based  on  random 
noise  can  only  qualitatively  explain  the  noise  effect 
on  the  experimental  data. 

Figure  10(b)  shows  the  experimental  results  of  re¬ 
constructed  |ia(max)  versus  the  total  number  of 
source— detector  pairs  obtained  from  3-D  imaging 
(circles)  and  ultrasound-assisted  imaging  (stars). 
The  upper  curve  is  the  fitting  result  of  the  circles,  and 
the  lower  curve  is  the  result  of  the  stars.  In  both 
cases,  the  reconstructed  values  were  decreased  when 
the  total  number  of  pairs  was  reduced.  However, 
the  reconstructed  values  were  more  accurate  when 
the  target  depth  information  was  available,  and  the 
improvement  on  average  was  15%.  The  improve¬ 
ment  was  more  dramatic  when  the  total  pairs  were 
less. 

For  sparse  arrays  with  total  source- detector 
pairs  less  than  140,  the  reconstructed  |io(maX)  could 
be  increased  if  the  iterations  were  significantly  in¬ 
creased.  However,  the  image  artifact  level  of  the 
target  layer  and  the  noise  level  of  the  nontarget 
layers  were  increased  too.  With  the  assistance  of 
simulations,  we  offer  the  following  explanations  to 
the  increased  image  artifact  and  noise  level  prob¬ 
lem.  In  simulations,  the  iteration  was  stopped 
when  the  TLS  error  between  the  measurement  and 
the  linear  approximation  reached  the  noise  floor, 
which  was 

E.t  <J  =  i  Ti(j)2  +  £  [2n(j)y\ 

.  •  j  .  . J  ... 

where  n(J)  was  the  generated  random  noise  with  the 
standard  deviation  proportional  to  a  certain  percent 
of  the  mean  of  the  forward  data  set  for  each  array 


configuration.  When  the  object  function  reached  the 
noise  floor,  the  gradient 

-2(£/sd-WAjOT(W) 

-2 (UK  -  ~  WA|xg)(A|xa) 

(A|x/  Ap.0  +  l)2 


can  be  approximated  as 


-2  (lV)r(W)  -2(AQT(AQ(A|xa) 
Ap.0TA|xa  +  1  (Ap.aTAp,a  +  l)2  ’ 


where  N  is  the  noise  vector.  Therefore  the  search 
procedure  of  Eq.  (11)  is  more  random  and  noisy.  Be¬ 
cause  A|x „  =  A|xaM((j  +  pVg,  the  Ap.  updating  is 
more  ranuom  and  noisy.  (3  is  proportional  to  the 
square  of  the  gradient.  Continuous  iteration  when 
the  object  function  has  reached  the  noise  floor  may 
destroy  the  convergence  of  the  reconstruction.  We 
found  that  when  the  SNR  of  the  data  is  high,  for 
example,  high-contrast  cases,  continuous  iteration  in 
general  increases  reconstructed  p.a  and  sidelobes. 
However,  when  the  SNR  of  the  data  is  low ,  for  exam¬ 
ple,  low-contrast  cases,  continuous  iteration  does  not 
increase  the  reconstructed  |xa,  but  destroys  the  con- 
vfirfence  of  the  reconstruction  [see  Fig.  12(c)  below]. 


F.  Experimental  Results  of  Imaging  Parameters  Versus 
Total  Number  of  Source-Detector  Pairs 
The  imaging  parameters  measured  from  different  ar¬ 
ray  configurations  are  listed  in  Table  1.  Listed  first 
are  the  measured  parameters  at  target  layer  3. 
These  parameters  are  a  —  6-dB  beam  width  of  the 
image  lobe,  the  peak  sidelobe  level,  |ia(max),  and  the 
distance  in  the  x-y  plane  between  the  location  of 
A'a(m.ax)  aad  the  true  tarSet  location.  Next  are  the 
same  parameters  measured  at  nontarget  layer  4. 
The  increase  in  beam  width  was  negligible  for  the 
arrays  with  more  than  140  source— detector  pairs  and 
was  100%  broader  for  the  sparse  arrays  with  total 
pairs  less  than  this  number.  The  sidelobe  level  was 
progressively  increased  when  the  total  pairs  were 
reduced.  At  nontarget  layer  4,  when  the  total  pairs 
were  reduced  to  less  than  140,  the  measured  image 
lobes  were  broad  and  no  sidelobes  were  seen.  The 
agreement  between  the  measured  |ia(max)  location 
and  the  true  target  location  is  good  for  all  the  array 
configurations,  which  suggests  that  this  parameter  is 
not  sensitive  to  the  total  number  of  source- detector 
pairs  in  high-contrast  cases.  Table  1  next  lists  the 
strength  of  the  target  measured  at  nontarget  layer  2 
and  its  spatial  location.  Because  the  SNR  of  the 
sparse  array  measurement  was  lower,  the  target  ap¬ 
peared  at  layer  2  when  total  pairs  were  less  than  180. 
However,  in  all  cases,  the  target  mass  measured  at 
this  layer  was  much  smaller  than  that  obtained  at  the 
target  layer  and  nontarget  layer  4.  Finally,  Table  1 
lists  the  measured  imaging  parameters  when  the  tar¬ 
get  depth  was  available  to  optical  reconstruction. 
Compared  with  parameters  obtained  from  optical  im- 
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Fig.  11.  Experimental  images  of  |ia  reconstructed  from  a  total  of 
28  X  13  =  364  source- detector  pairs.  The  target  (p.a  =  0.10  cm" L) 
was' located  at-  (jc  =  0,-  y :=  0/2'  =’:2.5;  cm)  inside  the  Intralipid 
background.^  i(a) 'Reconstructed  at  target  layer  3.  The  mea- 
sur.ed  p.a(max7  was  0.063  cm"1  at  510  iterations,  and  its  location  was 
(?l  —  p.O,  y.  Edge  artifacts  were  observed  and  the  peak 

level  was  ji7  ;dB  from, the  ,|xa(raa^/  (b)  Reconstructed  |la  at  non- 
target  layer  'far,  ;The^measured  jia&ax)  was  0.0871  cm  1  at  510 
iterations,  and (its  location  was  {x  =,0 ,;y  =  0).  (c)  Reconstructed 

^  ‘at  the  r  target ;la^er  when’  only  tlie  target  depth  information 
Provided  by /ultraiound  was  used'.// £a(max)  =  0.107  cm"1  at  56 

ikati 

agiri'g1  orilK^Ke^6-dB  beam  width  was  improved  by 
24%  on  average^  and  the  reconstruction  speed  was 
approximately;  10  times  faster;  however,  the  sidelobe 
was  3  dB  . worse* 

5.  Results  of  a  Lower-Contrast  Target  Case 

X;i;jxperirT|^^^|;of^  Dense  Array 

To  study  the  effects  of  target  contrast  on  the  quality 

of:rthe  reconstructed  linage  for  each  array  configura- 
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Intralipid^p^pP^teT;  of  the  target  was  located  at 
(jcT==  cm).  ■  .The  centers  of  the 

imaging  voxels.'. ifef^ere.  .0,5, -1.5, .  2.5,  and  3.5  cm  for 
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Fig.  12.  Experimental  images  of  |ia  at  target  layer  3  recon¬ 
structed  from  a  total  of  24  X  5  =  120  source-detector  pairs,  (a) 
Reconstructed  (la  at  target  layer  3  with  510  iterations.  The  mea¬ 
sured  |io(max)  was  0.049  cm"1,  which  was  49%  of  the  true  target  pa, 
and  its  location  was  (*  =  0.5,  y  =  - 1.0),  which  was  displaced  from 
the  true  target  location  by  1.11  cm.  Image  artifacts  were  ob¬ 
served,  and  the  peak  was  —3  dB  from  the  |la(mnx)-  (b)  —6-dB 
contour  plot  of  (a),  (c)  Reconstructed  |i0  at  target  layer  3  with 

1500  iterations.  The  peak  artifact  was  5  dB  higher  than  the 
image  lobe,  (d)  Reconstructed  |ia  at  target  layer  3  (56  iterations) 
with  only  a  priori  target  depth  information  provided  by  ultra¬ 
sound.  The  measured  (ia(maX|  was  0.074  cm  \  and  its  location 
was  ( x  =  0.5,  y  =  -0.5).  Image  artifacts  were  observed,  and  the 
peak  was  5  dB  from  the 


layers  1,  2,  3,  and  4,  respectively.  Figure  11  shows 
the  images  of  the  lower-contrast  target  obtained  at 
target  layer  3  [Fig.  11(a)]  and  deeper  nontarget  layer 
4  [Fig.  11(b)].  The  images  were  reconstructed  from 
the  measurements  made  with  the  dense  array.  At 
target  layer  3,  the  measured  |ia(max)  was  0.063  cm  , 
which  was  approximately  63%  of  the  target  |ia.  The 
measured  spatial  location  of  |ia(max) was  (*  =  = 

-0.5),  which  agreed  reasonably  well  with  the  true 
target  location.  The  measured  -6-dB  beam  width 
was  1.83  cm,  which  was  144%  times  broader  than  the 
beam  width  of  the  high-contrast  case.  Edge  arti¬ 
facts  were  observed,  and  the  peak  level  was  —7  dB 
from  the  jLo(max).  At  nontarget  layer  4,  the  measured 
ha(max) was  0.0871  cm-1,  which  was  even  higher  than 
that  measured  at  the  target  layer.  Because  the  SNR 
of  the  data  was  lower  than  the  high-contrast  case,  the 
target  depth  estimate  was  poorer.  A  target  of 
ha  (max)  =  0.0543  cm'1  located  at  (x  =  0,  y  =  -0.5) 
was  observed  at  nontarget  layer  2,  and  its  mass  was 
much  smaller  than  that  obtained  at  the  target  layer 
and  nontarget  layer  4. 

Figure  11(c)  is  the  reconstructed  |ia  at  the  target 
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layer  fr&m  only  the  dense  array  measurements.  A 
total  of  56  iterations  were  used  to  obtain  (io(max)  = 
0.107  cm  \  and  the  reconstruction  was  approxi¬ 
mately  ten  times  faster  than  that  without  the  depth 
information.  The  spatial  resolution  was  8%  better 
than  that  obtained  from  Fig.  11(a),  and  the  -6-dB 
beam  width  was  1.67  cm.  The  contrast  resolution 
was  2  dB  worse. 

B.  Experimental  Results  of  a  Sparse  Array 

The  imaging  quality  of  sparse  arrays  decreased. 
Figure  12(a)  is  an  image  of  the  same  target  (fxa  =  0.10 
cm  ’)  reconstructed  from  measurements  made  with 
the  24  x  5  sparse  array.  The  measured  |ia(max)  was 
0.049  cm  \  which  was  49%  of  the  true  value.  The 
-6-dB  contour  plot  is  shown  in  Fig.  12(b).  The  mea¬ 
sured  spatial  location  of  |ia(max)  was  (x  =  0.5,  y  = 
—  1.0),  which  was  displaced  from  the  true  target  lo¬ 
cation  by  1.11  cm  in  radius.  The  measured  -6-dB 
beam  width  was  2.98  cm,  which  was  163%  broader 
than  that  measured  from  the  dense  array.  Side- 
lobes  were  abundant,  and  the  peak  value  was  -3  dB 
below  the  peak  of  the  image  lobe.  These  sidelobes 
would  produce  false  targets  in  the  image  if  no  a  priori 
information  about  the  target  locations  were  given. 

In  this  case,  continuous  iteration  did  not  increase 
the  target  strength  but  increased  the  sidelobe 
strength.  Figure  12(a)  was  obtained  at  510  itera¬ 
tions,  whereas  Fig.  12(c)  was  obtained  at  1500  itera¬ 
tions.  After  approximately  three  times  more 
iterations,  the  peak  of  the  artifact  was  5  dB  higher 
than  the  peak  of  the  image  lobe.  The  target  depth 
estimated  from  3-D  images  was  worse  at  1500  itera¬ 
tions  than  that  at  the  510  iterations.  The  measured 
target  strengths  at  nontarget  layer  2  were  0.064  and 
0.1661  cm-1  at  510  and  1500  iterations,  respectively, 
and  the  strengths  at  nontarget  layer  4  were  0.0597 
and  0.1087  cm-1,  respectively.  In  addition,  the 
background  noise  fluctuation  or  standard  deviation 
measured  at  nontarget  layer  1  was  increased  with  the 
iterations.  The  mean  and  the  standard  deviation  at 
510  iterations  were  0.024  and  0.0041  cm-1,  whereas 
these  values  at  1500  iterations  were  0.0267  and 
0.0104  cm-1. 

As  shown  in  Fig.  12(d),  the  reconstructed  image 
improved  a  lot  when  the  target  depth  was  given. 
The  maximum  strength  was  0.074  cm-1,  and  its  lo¬ 
cation  was  (0.5,  —0.5).  The  —6-dB  beam  width  was 
2.63  cm,  which  was  12%  better  than  that  obtained 
from  Fig.  12(a).  The  sidelobe  was  -6  dB  from  the 
peak,  which  was  improved  by  3  dB  compared  with 
Fig.  12(a). 

C.  Experimental  Results  of  Reconstructed  (ia(max)  versus 
Total  Number  of  Source-Detector  Pairs 

With  the  iteration  number  normalized  to  the  dense 
array  case,  we  measured  target  strengths  at  target 
layer  3  and  nontarget  layer  4  for  all  sparse  array 
configurations.  Figure  13  shows  the  experimental 
data  points  (circles)  of  measured  |ia(max)  versus  the 
total  number  of  source- detector  pairs  obtained  at 
target  layer  3  [Fig.  13(a)]  and  nontarget  layer  4  [Fig. 


50  100  150  200  250  300  350  400 

total  number  of  soiree-detector  pairs 


(a) 


total  source-detector  pairs 

(b) 

Fig.  13.  Low-contrast  target  case,  (a)  Reconstructed  |io(max)  ver¬ 
sus  total  source- detector  pairs  with  the  target  depth  available 
(stars)  and  the  curve-fitting  results  (upper  curve).  Reconstructed 
M-a(max)  versus  total  source- detector  pairs  measured  at  target  layer 
(circles)  and  the  curve-fitting  results  (lower  curve)  and  (b)  at 
deeper  nontarget  layer  4. 


13(b)].  The  two  curves  were  the  fitting  results  of 
experimental  points  when  we  used  second-order  poly¬ 
nomials.  In  both  layers,  the  (la(max)  values  were  de¬ 
creased  when  the  total  number  of  source- detector 
pairs  was  reduced.  The  reconstructed  target 
strengths  were  reduced  to  less  than  60%  when  the 
total  pairs  were  less  than  156  and  140  for  target  layer 
3  and  nontarget  layer  4,  respectively.  The 
ultrasound-assisted  reconstruction  results  are  shown 
in  Fig.  13(a)  (stars),  and  the  average  reconstructed 
Fa(max)  f°r  all  array  configurations  was  0.085  cm-1. 
Compared  with  the  average  of  0.055  cm-1  obtained 
from  optical  imaging  only  at  the  target  layer,  a  30% 
improvement  was  achieved. 

D.  Experimental  Results  of  Imaging  Parameters  versus 
Total  Number  of  Source-Detector  Pairs 

The  measured  imaging  parameters  obtained  from  dif¬ 
ferent  array  configurations  are  listed  in  Table  2. 
Similar  to  Table  1,  listed  are  the  measured  parame¬ 
ters  at  target  layer  3  and  nontarget  layer  4.  The 
increase  in  beam  width  at  the  target  layer  was  neg¬ 
ligible  for  the  arrays  with  more  than  156  source- 
detector  pairs  and  was  more  than  50%  for  the  sparse 
arrays  with  total  pairs  less  than  this  number.  The 
sidelobes  were  progressively  worse  when  the  total 
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pairs  were  reduced.  At  the  target  layer,  the  distance 
between  the  |ia(max)  location  and  the  true  target  loca¬ 
tion  was  more  than  1  cm  for  the  arrays  with  total 
pairs  less  than  140.  Table  2  then  lists  the  measured 
peak  image  lobe  at  nontarget  layer  2  and  its  location. 
The  target  appeared  at  nontarget  layer  2  in  all  array 
configurations  because  of  the  lower  SNR  of  the  data. 
However,  the  target  mass  observed  at  this  layer  was 
much  smaller  than  that  at  layers  3  and  4  in  all  cases. 
The  last  row  in  Table  2  shows  the  measured  imaging 
parameters  when  the  target  depth  was  available  to 
optical  reconstruction.  Compared  with  parameters 
obtained  from  optical  imaging  only,  the  -6-dB  beam 
width  was  improved  by  41%  on  average,  and  the  re¬ 
construction  speed  was  approximately  ten  times  fast¬ 
er;  however,  the  sidelobe  was  1  dB  worse. 

6.  Discussion 

In  addition  to  the  total  number  of  source- detector 
pairs,  the  measured  imaging  parameters  are  also  re¬ 
lated  to  other  system  parameters,  for  example,  mod¬ 
ulation  frequency  and  system  noise.  The  140-MHz 
modulation  frequency  chosen  in  this  study  is  a  typical 
frequency  used  by  many  research  groups.  Our  sys¬ 
tem  noise  level,  including  both  coherent  and  incoher¬ 
ent,  is  less  than  10  mV  peak  to  peak,  which  is 
sufficiently  low.  Therefore  the  results  we  obtained 
are  pertinent  to  3-D  imaging  using  similar  system 
parameters  and  reflection  geometry. 

In  this  study  the  target  absorption  coefficient  was 
reconstructed  from  measurements.  Similar  studies 
can  be  done  for  the  scattering  coefficient  as  well.  To 
reconstruct  the  scattering  coefficient,  we  can  use  the 
scattering  weight  matrix  derived  by  O’Leary11  to  re¬ 
late  the  medium  scattering  variations  with  the  mea¬ 
surements.  Simultaneous  reconstruction  of  both 
absorption  and  scattering  coefficients  in  reflection  ge¬ 
ometry  is  also  possible,  provided  that  the  absorption 
and  scattering  weight  matrices  are  regulated  care¬ 
fully.  Because  the  eigenvalues  of  the  two  matrices 
are  significantly  different,  good  regulation  schemes 
are  needed  to  balance  the  reconstructed  absorption 
and  scattering  coefficients  at  each  iteration.  This 
subject  is  one  of  our  topics  for  further  study. 

In  this  study  the  ultrasound-assisted  optical  recon¬ 
struction  was  demonstrated  at  a  particular  target 
layer.  Similar  studies  can  be  done  with  multiple 
targets  located  at  different  layers.  In  the  multiple 
target  case,  we  can  attribute  measured  perturbations 
to  more  layers  instead  of  a  single  layer.  However, 
the  improvements  in  reconstructed  optical  properties 
and  reconstruction  speed  may  be  less  than  that  of  the 
single-layer  case. 

Ultrasound  has  good  imaging  capability,  and  it  can 
detect  small  lesions  of  a  few  millimeters  in  size. 
However,  its  specificity  in  cancer  detection  is  not  high 
as  a  result  of  overlapping  characteristics  of  benign 
and  malignant  lesions.  NIR  imaging  has  high  spec¬ 
ificity  in  cancer  detection;  however,  it  suffers  low  res¬ 
olution  and  lesion  location  uncertainty  because  of  the 
diffused  nature  of  the  NIR  light.  The  hybrid  imag¬ 
ing  that  combines  ultrasound  imaging  capability  and 


NIR  contrast  has  a  great  potential  to  overcome  defi¬ 
ciencies  of  either  method.  As  we  reported  in  the 
paper,  the  target  depth  information  can  significantly 
improve  the  accuracy  of  the  reconstructed  optical  ab¬ 
sorption  coefficient  and  reconstruction  speed.  In  ad¬ 
dition  to  use  of  a  priori  target  depth  information,  the 
target  spatial  distribution  provided  by  ultrasound 
can  be  used  in  optical  reconstruction  as  well.16 
With  the  localized  spatial  and  temporal  target  in¬ 
formation,  the  accuracy  of  the  reconstructed  optical 
properties  and  the  reconstruction  speed  can  be  im¬ 
proved  further.  To  demonstrate  this,  we  need  an 
ultrasound  imaging  transducer  located  at  the  mid¬ 
dle  portion  of  the  probe.  We  are  currently  pursu¬ 
ing  this  study.23-24 

In  this  study  the  targets  of  different  contrasts  were 
located  at  the  center  position.  We  have  also  done 
studies  with  targets  of  different  contrasts  located  at 
off-center  positions.  For  an  off-center  target  case, 
the  effective  number  of  source- detector  pairs  is  less 
than  that  of  the  on-center  target  case  because  mea¬ 
surements  from  certain  source- detector  pairs  do  not 
contain  much  information  about  the  target.  For  ex¬ 
ample,  if  a  target  is  placed  at  (x  =  2,  y  -  2,  z  =  3.0 
cm),  the  measurements  of  source- detector  pairs  at 
the  opposite  corner  of  the  probe  contribute  little  to  the 
image  reconstruction.  In  one  study,  targets  of  high 
and  low  contrast  were  located  at  (x  =  2,  y  =  2,  z  =  3.0 
cm).  The  reconstructed  maximum  absorption  coef¬ 
ficient  at  the  target  layer  was  related  more  to  the 
total  neighbor  source- detector  pairs.  However,  the 
maximum  value  at  a  deeper  nontarget  layer  was  re¬ 
lated  more  to  the  total  source-detector  pairs  and  was 
decreased  with  the  reduction  of  total  pairs.  Because 
the  photons  originated  from  distant  sources  and  de¬ 
tected  by  distant  detectors  experience  longer  and 
more  diffused  scattering  paths,  they  are  likely  to  in¬ 
teract  with  the  off-center  target  and  contribute  to  the 
absorption  estimate  at  the  deeper  layer.  In  the 
same  study,  the  measured  sidelobe  levels  were  4.5 
and  2.0  dB  poorer  on  average  compared  with  the 
on-center  high-  and  low-target  cases,  respectively. 
The  beam  widths  were  comparable  to  those  measured 
from  on-center  cases. 

7.  Summary 

The  relationship  between  the  total  number  of  source- 
detector  pairs  and  the  imaging  parameters  of  a  re¬ 
constructed  absorption  coefficient  was  evaluated 
experimentally.  A  frequency-domain  system  of  a 
140-MHz  modulation  frequency  was  used  in  the  ex¬ 
periments.  Reconstruction  at  a  selected  target 
depth  with  a  priori  depth  information  provided  by 
ultrasound  was  demonstrated.  The  results  have 
shown  that  the  reconstructed  absorption  coefficient 
and  the  spatial  resolution  of  the  absorption  image 
were  decreased  when  the  total  number  of  source- 
detector  pairs  was  reduced.  More  than  160  source- 
detector  pairs  were  needed  to  reconstruct  the 
absorption  coefficient  within  60%  of  the  true  value 
and  spatial  resolution  comparable  to  that  obtained 
with  the  dense  array.  The  contrast  resolution  was 
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poorer  in  general  because  of  edge  artifacts  and  could 
be  worse  if  significant  larger  iteration  numbers  are 
used  for  reconstruction.  The  error  in  target  depth 
estimated  from  3-D  optical  images  was  approxi¬ 
mately  1  cm.  With  the  a  priori  target  depth  infor¬ 
mation  provided  by  ultrasound,  the  reconstruction 
can  be  done  at  a  selected  depth.  Because  the  un¬ 
knowns  were  reduced  significantly,  the  reconstruc¬ 
tion  speed  was  approximately  ten  times  faster  than 
that  without  depth  information.  In  addition,  the  ac¬ 
curacy  of  the  reconstructed  absorption  coefficient  was 
improved  by  15%  and  30%  on  average  for  high-  and 
low-contrast  cases,  respectively.  Furthermore,  the 
measured  -6-dB  beam  width  was  improved  by  24% 
and  41%  for  high-  and  low-contrast  cases,  respec¬ 
tively.  The  sidelobe  was  3  and  1  dB  poorer  for  high- 
and  low-contrast  cases  because  the  measurement 
noise  was  lumped  to  single-layer  reconstruction  in¬ 
stead  of  multiple  layers. 

In  conclusion,  ultrasound-assisted  3-D  optical  im¬ 
aging  has  shown  promising  results  to  overcome  the 
problems  associated  with  the  reconstruction  by  use  of 
diffusive  waves.  With  the  target  depth  information 
provided  by  ultrasound,  the  reconstructed  absorption 
coefficient  was  more  accurate  and  the  reconstruction 
speed  was  much  faster. 
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Abstract 

We  have  constructed  a  near  real-time  combined  imager  suitable  for  simultaneous  ultrasound  and 
near  infrared  (NIR)  diffusive  light  imaging  and  co-registration.  The  imager  consists  of  a 
combined  hand-held  probe  and  the  associated  electronics  for  data  acquisition.  A  two- 
dimensional  ultrasound  array  is  deployed  at  the  center  of  the  combined  probe  and  12  dual 
wavelength  laser  source  fibers  (780  nm  and  830  nm)  and  8  optical  detector  fibers  are  deployed  at 
the  periphery.  We  have  experimentally  evaluated  the  effects  of  missing  optical  sources  in  the 
middle  of  the  combined  probe  upon  the  accuracy  of  the  reconstructed  optical  absorption 
coefficient,  and  assessed  the  improvements  of  reconstructed  absorption  coefficient  with  the 
guidance  of  the  co-registered  ultrasound.  The  results  have  shown  that  when  the  central 
ultrasound  array  area  is  in  the  neighborhood  of  2x2  cm2,  which  corresponds  to  the  size  of  most 
commercial  ultrasound  transducers,  the  optical  imaging  is  not  affected.  The  results  have  also 
shown  that  the  iterative  inversion  algorithm  converges  very  fast  with  the  guidance  of  a  priori 
three-dimensional  target  distribution,  and  only  one  iteration  is  needed  to  reconstruct  accurate 
optical  absorption  coefficient. 
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1.  Introduction 

Ultrasound  is  used  extensively  for  differentiation  of  cysts  from  solid  lesions  in  breast 
examinations  and  it  is  routinely  used  in  conjunction  with  mammography.  Ultrasound  can  detect 
breast  lesions  a  few  mm  in  size.1  However,  its  specificity  in  breast  cancer  detection  is  not  high 
as  a  result  of  overlapping  characteristics  of  benign  and  malignant  lesions.  Optical  imaging 
based  on  diffusive  near-  infrared  (NIR)  light  has  the  great  potential  to  differentiate  tumors  from 
normal  breast  tissues  through  determination  of  tissue  parameters,  such  as  blood  volume,  blood 
02  saturation,  tissue  light  scattering,  water  concentration,  and  the  concentration  and  lifetime  of 
exogenous  contrast  agents.4"12  As  a  potential  diagnostic  tool,  however,  NIR  diffusive  light 
imaging  suffers  from  low  spatial  resolution  and  lesion  location  uncertainties  due  to  intense  light 
scattering  in  tissue. 

Most  NIR  imaging  reconstruction  algorithms  are  based  on  tomographic  inversion  techniques.13"19 
Reconstruction  of  tissue  optical  properties  in  general  is  underdetermined  and  ill-posed  because 
the  total  number  of  unknown  optical  properties  always  exceeds  the  number  of  measurements  and 
the  perturbations  produced  by  the  heterogeneities  are  much  smaller  than  the  background  signals. 
In  addition,  the  inversion  reconstruction  algorithms  are  very  sensitive  to  measurement  noise  and 
model  errors. 

Our  group  and  others  have  introduced  a  novel  hybrid  imaging  method  that  combines  the 
complementary  features  of  ultrasound  and  near-infrared  diffusive  light  imaging.20'24  The  hybrid 
imaging  obtains  co-registered  ultrasound  and  NIR  diffusive  light  images  through  simultaneous 
deployment  of  an  ultrasound  array  and  NIR  source  detector  fibers  on  the  same  probe.  ’  ’  Co- 
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registration  permits  joint  evaluation  of  acoustic  and  optical  properties  of  breast  lesions,  and 
enables  the  use  of  lesion  morphology  provided  by  high-resolution  ultrasound  to  improve  the 
lesion  optical  property  estimate.  With  the  a  priori  knowledge  of  lesion  location  and  shape 
provided  by  co-registered  ultrasound,  NIR  imaging  reconstruction  can  be  localized  within 
specified  3-D  regions.  As  a  result,  the  reconstruction  is  over-determined  because  the  total 
number  of  unknown  optical  properties  is  reduced  significantly.  In  addition,  the  reconstruction  is 
less  sensitive  to  noise  because  the  convergence  can  be  achieved  with  a  small  number  of 
iterations. 

The  clinical  use  of  the  combined  diagnosis  relies  on  the  co-registration  of  both  ultrasound  and 
NIR  sensors  at  the  probe  level.  Conventional  ultrasound  pulse  echo  imaging  requires  an  imaging 
transducer  be  located  on  top  of  the  target,  while  NIR  diffusive  light  imaging  is  feasible  when  the 
optical  source  and  detector  fibers  are  distributed  at  the  periphery  of  the  ultrasound  transducer. 
However,  the  effects  of  missing  optical  sources  in  the  middle  of  the  combined  probe  upon  the 
accuracy  of  the  reconstructed  optical  properties  have  to  be  evaluated.  In  addition,  the 
improvements  of  reconstructed  optical  properties  with  the  guidance  of  the  co-registered 
ultrasound  need  to  be  quantitatively  assessed.  Furthermore,  real-time  data  acquisition  is 
necessary  to  avoid  errors  in  co-registration  caused  by  patient  motion  during  the  clinical 
experiments.  In  this  paper,  we  report  our  experimental  results  on  optimal  probe  configuration 
and  we  quantify  the  improvements  on  reconstructed  optical  properties  using  a  combined  probe. 
We  also  demonstrate  simultaneous  combined  imaging  with  a  near  real-time  imager. 


2.  NIR  diffusive  wave  imaging 


4 


We  have  used  Bom  approximation  to  relate  the  scattered  field  Usc(r,a>)  measured  at  the  probe 
surface  to  absorption  variations  in  each  volume  element  within  the  sample.  In  the  Born 
approximation,  the  scattered  wave  originated  from  a  source  at  r„  and  measured  at  vdi  can  be 

related  to  the  medium  absorption  heterogeneity  A jia  (rv> )  at  rvj  by 

=  }/V*l 

where  M  is  the  total  number  of  source-detector  pairs,  N  is  the  total  number  of  imaging  voxels, 
and  w0  =  G(rvj ,  rdi ,  Q))Uinc (rvj , r„. ,  (0)vA r* •/  D  is  the  weight  matrix  given  in  Ref.  18. 
G(rVj,rid,0))  and  Uinc (rvy. , rt/ , co) are  Green  function  and  incident  wave,  respectively.  (0  is  the 

modulation  frequency  and  Dis  the  average  or  background  diffusion  coefficient,  which  is  the 
average  value  over  the  background  or  whole  tissue. 

With  M  measurements  obtained  from  all  possible  source-detector  pairs  in  the  planar  array,  we 

can  solve  N  unknowns  of  fiu  by  inverting  the  above  matrix  equation.  In  general,  the 
perturbation  Eq.  (1)  is  underdetermined  (M<N)  and  ill-posed. 

NIR  imaging  by  itself  generally  has  poor  depth  discrimination.  However,  ultrasound  is  excellent 
in  providing  accurate  target  depth.  Once  the  target  depth  is  available  from  co-registered 
ultrasound,  we  can  set  Ajua  of  non-target  depth  equal  to  zero.  This  implies  that  all  the  measured 
perturbations  originate  from  the  particular  depth  that  contains  the  target.  Since  the  number  of 
unknowns  is  reduced  significantly,  the  reconstruction  converges  very  fast.  In  Ref  [23],  we  have 
reported  that  with  a  priori  target  depth  provided  by  ultrasound,  the  accuracy  of  the  reconstructed 
fla  has  been  improved  by  15%-30%  on  average  and  the  speed  of  reconstruction  has  been 
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improved  by  an  order  of  magnitude.  In  this  paper,  we  furthermore  demonstrate  that  with  the 
three-dimensional  target  distribution  provided  by  co-registered  ultrasound,  the  accuracy  of 
reconstructed  //„  and  the  reconstruction  speed  can  be  further  improved. 


To  solve  the  unknown  optical  properties  of  Eq.  (1),  we  have  used  the  total  least-squares  (TLS) 
method  25-26  to  iteratively  invert  the  equation.  The  TLS  performs  better  than  other  least-square 
methods  when  the  measurement  data  are  subject  to  noise  and  the  linear  operator  W  contains 
errors.  We  have  found  that  the  TLS  method  provides  more  accurate  reconstructed  optical 
properties  than  other  least  square  methods  and  we  have  adopted  TLS  in  solving  inverse 
problems.  It  has  been  shown  in  Ref.  [27]  that  the  TLS  minimization  is  equivalent  to  the 
following  minimization  problem. 


min 


(2) 


where  X  represents  unknown  optical  properties.  The  conjugate  gradient  technique  has  been 


employed  to  iteratively  solve  Eq.  (2). 


3.  Methods 

A.  Combined  probe  and  imaging  geometry 

There  are  four  basic  requirements  to  guide  the  design  of  the  combined  probe.  First,  reflection 
geometry  is  preferred  because  conventional  ultrasound  scan  is  performed  with  this  geometry. 
Second,  ultrasound  array  needs  to  occupy  the  center  of  the  combined  probe  for  coherent  imaging. 
Third,  NIR  sources  and  detectors  have  to  be  distributed  at  the  periphery.  Since  photon 
propagation  distribution  exhibits  a  “banana”  shape,  imaging  of  the  tissue  volume  underneath  the 
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probe  is  feasible  even  through  there  are  no  sources  and  detectors  deployed  in  the  central  portion 
of  the  probe.  Fourth,  the  minimum  source-detector  separation  should  be  larger  than  1  cm  for 
diffusion  approximation  to  be  valid,  and  the  maximum  separation  should  be  ~  8-9  cm  for 
effectively  probing  depths  of  3-4  cm. 

Based  on  these  requirements  we  have  deployed  12  dual-wavelength  optical  source  fibers  and  8 
detector  fibers  over  a  9x9  cm2  probe  area  (see  Figure  1).  The  minimum  and  maximum  source- 
detector  separations  in  the  configuration  are  1.4  cm  and  8  cm,  respectively.  To  study  the  effect 
of  the  central  “optical  hole”  upon  the  accuracy  of  the  reconstructed  optical  properties,  we  have 
compared  the  reconstruction  results  with  an  extra  center  source  and  without  the  center  source. 
The  configuration  without  the  center  source  corresponds  to  a  2  x  2  cm  hole  area.  We  further 
moved  the  non-center  12  sources  and  8  detectors  toward  periphery  by  leaving  a  3x3  cm  hole 
area  in  the  middle.  Figure  2  shows  the  picture  of  a  combined  probe  with  the  3x3  cm  central 
area  occupied  by  an  ultrasound  array.  The  ultrasound  array  consists  of  64  elements  made  of  1.5 
mm  diameter  piezoelectric  (PZT)  transducers  (Valpey  Fisher  Inc).  The  transducers  are  deployed 
in  a  rectangular  matrix  with  4  mm  spacing  in  both  x  and  y  directions.  The  center  frequency  of 
the  transducer  is  6  MHz  and  the  bandwidth  is  40%.  The  transducers  are  made  from  the  same 
piece  of  PZT  material.  Therefore,  the  gain  difference  among  different  transducers  is  less  than  3 
dB.  12  dual-wavelength  optical  laser  diode  sources  (760  nm  and  830  nm)  and  8  Photo  Multiplier 
(PMT)  detectors  are  coupled  to  the  probe  through  optical  fibers,  which  are  deployed  at  the 
periphery  of  the  2-D  ultrasound  array.  This  hybrid  array  deployment  compromises  ultrasound 
coherent  imaging  and  NIR  diffusive  light  imaging  characteristics. 
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The  9x9x4  cm3  image  volume  underneath  the  probe  is  discretized  into  voxels  of  size  0.4 x 0.4 
x  1  cm3.  There  is  a  trade-off  between  the  accurate  estimation  of  the  weight  matrix  W  and  the 
voxel  size.  Since  Wt]  is  a  discrete  approximation  of  the  integral 


v 


the  total  number  of  reconstructed  unknowns  will  increase  dramatically  with  the  decreasing  voxel 
size.  Furthermore,  the  rank  of  W  does  not  increase  in  the  same  order  as  the  total  number  of 
voxels  when  the  voxel  size  decreases.  This  suggests  that  neighboring  Wtj  s  are  correlated  when 

the  voxel  size  is  smaller  and  further  decrease  in  voxel  size  will  not  add  more  independent 
information  to  the  weight  matrix.  We  have  found  that  0.4 x  0.4x1  cm3  voxel  size  is  a  good 
compromise.  Therefore,  we  have  used  this  voxel  size  in  image  reconstructions  reported  in  this 
paper. 


B.  Experimental  Systems 
NIR  imaging  system 

We  have  constructed  a  NIR  frequency  domain  imaging  system.  The  block  diagram  of  the  system 
is  shown  in  Figure  3.  This  system  has  12  dual -wavelength  source  channels  and  8  parallel 
receiving  channels.  On  the  transmission  part,  12  pairs  of  dual  wavelength  (780nm  and  830nm) 
laser  diodes  are  used  as  light  sources,  and  their  outputs  are  amplitude  modulated  at  140.000 
MHz.  Each  one  of  the  12  optical  combiners  (OZ  Optics  Inc.)  looks  like  a  Y  adapter,  guiding  the 
emission  of  two  diodes  of  different  wavelengths  through  the  same  thin  optical  fiber  (about  0.2 
mm  in  diameter).  To  reduce  noise  and  interference,  an  individual  driving  circuit  is  built  for  each 
diode.  As  a  laser  diode  works  in  series,  a  control  board  that  interprets  instructions  from  a  PC  is 
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used  to  coordinate  operations  of  associated  components.  When  a  single  transmission  channel  is 
selected,  it  turns  on  corresponding  driving  circuit  so  that  a  DC  driving  current  can  be  set  up  for 
the  diode.  At  the  same  time,  a  selected  signal  is  sent  to  a  RF  switching  unit,  which  distributes  a 
RF  signal  to  the  right  channel  to  modulate  the  optical  output.  On  the  reception  part,  8  Photo 
Multiplier  Tubes  (PMTs)  are  employed  to  detect  diffusely  reflected  light  from  turbid  media. 
Each  PMT  is  housed  in  a  sealed  aluminum  box,  shielding  both  environment  lights  and 
electromagnetic  fields,  while  an  optical  fiber  (3  mm  in  diameter)  couples  NIR  light  from  the 
detection  point  to  the  reception  window  of  the  PMT.  The  electrical  signal  converted  from  the 
optical  input  is  generally  very  weak  and  rather  high  in  frequency,  so  high  gain  amplification  and 
frequency  transform  are  necessary  before  it  can  be  sampled  by  an  A/D  board  inside  the  PC.  We 
have  built  8  parallel  heterodyne  amplification  channels  in  order  to  measure  the  response  of  all 
detectors  simultaneously,  to  thus  reduce  the  data  acquisition  time.  Each  amplification  channel 
consists  of  a  RF  amplifier  (40dB),  a  mixer  where  RF  signal  (OSC1,  140.000  MHz)  is  mixed  with 
a  local  oscillator  (OSC2,  140.020  MHz)  and,  a  band  pass  filter  centered  at  20  KHz,  and  a  low 
frequency  amplifier  of  30dB.  The  heterodyned  two-stage  amplification  scheme  helps  suppress 
wide  band  noises  efficiently.  A  reference  signal  of  20  KHz  is  also  generated  by  directly  mixing 
OSC1  and  OSC2,  which  is  necessary  for  retrieving  phase  shifts.  Eight  detection  signals  and  one 
reference  are  sampled,  converted,  and  acquired  into  the  PC  simultaneously,  in  which  the  Hilbert 
transform  is  used  to  compute  the  amplitude  and  phase  of  each  channel.  The  entire  data 
acquisition  takes  less  than  one  minute,  which  is  fast  enough  to  acquire  data  from  patients. 


One  of  the  challenges  encountered  in  designing  a  NIR  imaging  system  is  the  huge  dynamic  range 
of  signals  received  at  various  source-detector  distances.  For  example,  for  a  semi-infinite 
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phantom  made  of  0.5%  Intralipid  solution,  the  amplitude  measured  at  1  cm  away  from  a  source 
is  about  5000  times  larger  than  that  at  8  cm  separation.  In  addition,  the  perturbation  due  to  an 
embedded  heterogeneity  with  optical  properties  similar  to  a  tumor  is  normally  a  few  percent  of 
the  background  signal.  As  a  result,  a  reflection  mode  NIR  imaging  system  should  have  at  least 
120  dB  dynamic  range  to  probe  a  target  up  to  4  cm  in  depth.  It’s  very  hard  to  build  amplifiers 
working  linearly  over  such  a  wide  dynamic  range.  We  overcome  this  difficulty  by  implementing 
two-level  source  outputs.  The  DC  output  of  a  laser  diode  is  controlled  by  adjusting  its  feedback 
loop,  while  the  RF  signal  is  switched  simultaneously  via  a  two-step  attenuator  (no  attenuation  or 
30  dB  attenuation).  When  the  source  and  detector  are  close  to  each  other,  the  source  is 
controlled  to  yield  low-level  output.  When  the  separation  gets  larger,  30  dB  higher  output  level 
should  be  used.  With  this  two-level  source  scheme,  our  system  achieved  fairly  good  linearity 
over  a  wide  range  of  source-detector  separations  (from  1.5  to  8  cm). 


Since  parameters  of  individual  laser  diode  and/or  PMT  vary  considerably  from  one  to  another, 
we  have  to  calibrate  the  gain  and  phase  shift  for  each  channel.  A  set  of  measurements  obtained 
from  all  source-detector  pairs  placed  on  the  boundary  of  a  homogeneous  medium  are 

Aap  and  (j)ap,  a  =  1,2 -  1,2,..., n. 


Here,  amplitude  Aap  and  phase  (f)ap  are  related  to  source  a  and  detector  fi ,  while  m  and  n  are 

the  total  number  of  sources  and  detectors  respectively.  From  the  diffusion  theory,  we  can  obtain 
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the  following  set  of  equations. 


Aap  =  Is(a)Id(fi) 


exp  (-kip#) 

Pip 


l  <t>ap=<Ps(<X)  +  <PAfi)  +  krPaP 
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in  which  /,(«) and (ps(oc) are  relative  gain  and  phase  delay  associated  with  source  channel  a, 

Id(fi) and (pd(P)  are  similar  quantities  associated  with  detector  channel  /?.  pup  is  the 

corresponding  separation,  and  kr  +  jkt  is  the  complex  wave  number.  We  obtain  the  following 

set  of  linear  equations  by  taking  logarithm  of  the  above  equations  related  to  amplitude 

|  log  {plpAap )  =  log(/,  ( a ))  +  log  (Id  (/?))  -  ktpap  (3) 

|  !>ap  =  <Ps  («)  +  <Pd(P)  +  KP«p 

Although  the  optical  properties  of  the  calibration  medium  are  known  in  advance,  we  leave  the 
wave  number  as  unknown.  u  main  reason  that  we  leave  kr  and  k,  as  unknowns  is  that  we 


want  to  include  this  calibration  method  in  our  clinical  experiment  procedures.  That  means  in 
vivo  calibration  with  the  probe  on  a  human  breast.  In  this  case,  the  optical  properties  of  the 
background  medium  are  completely  unknown.  We  have  verified  our  calibration  method  by 
comparing  the  best  fitted  k ’s  with  real  values.  The  results  of  using  0.5~0.8%  Intralipid  solutions 
always  gave  scattering  and  absorption  coefficients  with  a  rather  good  accuracy.^  With  the  two 
unknown  wave  numbers  included,  the  total  number  of  unknowns  is  2(m  +  n  +  2),  which  is 
generally  far  smaller  than  the  number  of  measurements  m  x  n .  Consequently,  Eq.  (3)  is  over¬ 
determined.  We  can  solve  all  I s(oc) , Id(J3) ,<ps(a)  and  (PjiP )  terms  as  well  as  two  unknown 
wave  numbers  in  a  least  square  sense.  Then  all  measurements  can  be  calibrated  accordingly. 

The  results  of  amplitude  Aap  =  CXP(  kfap)- and  phase  <pap  =  krpap  after  calibration  are  shown 

Pap 

in  Figure  4.  As  one  can  see,  the  calibrated  amplitude  (log  pjA^and  phase  from  various 
source-detector  pairs  change  linearly  with  distance. 


Ultrasound  system 
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The  ultrasound  system  diagram  is  shown  in  Figure  5,  and  the  system  consists  of  64  parallel 
transmission  and  receiving  channels.  Each  transmission  circuit  can  generate  a  high  voltage  pulse 
of  200  ns  duration  (6  MHz)  with  125  volts  peak  to  peak  to  the  connected  transducer.  Each 
receiving  circuit  has  two-stage  amplifiers  followed  by  an  A/D  converter  with  40  MHz  sampling 
frequency.  The  amplifier  gain  can  be  controlled  based  on  the  target  strength.  A  group  of 
transmission  channels  can  be  addressed  simultaneously  to  transmit  pulses  from  neighbor 
transducers  with  specified  delays  and,  therefore,  to  focus  the  transmission  beam.  The  retuned 
signals  can  be  simultaneously  received  by  a  group  of  transducers  and  the  signals  can  be  summed 
with  specified  delays  to  form  a  receiving  beam. 

The  data  acquisition  procedure  is  the  following.  The  first  9-element  neighbor  subarray  (dashed 
rectangle  in  Figure  6)  from  the  64-element  transducer  array  and  the  corresponding  channels  are 
chosen,  and  then  the  transmission  delay  profiles  are  generated  in  the  computer  according  to  the 
pre-specified  focal  depth.  The  delay  profile  data  are  transferred  to  the  64-channel  delay  profile 
generator,  which  triggers  the  64  high  voltage  pulsers  as  well  as  the  receiving  channels.  The 
returned  ultrasound  signals  are  amplified  by  two  stage  amplifiers  and  sampled  by  A/D 
converters.  The  data  are  buffered  in  the  memories  and  are  read  by  the  computer  after  the  entire 
data  acquisition  is  completed.  The  second  subarray  (solid  rectangle  on  Figure  6)  is  chosen  and 
the  same  data  acquisition  process  is  repeated.  A  total  of  64  subarrays  is  used  in  the  data 
acquisition.  After  the  64-subarray  data  acquisition  is  completed,  the  data  stored  in  the  memories 
are  read  by  the  computer  for  image  formation.  The  entire  data  acquisition  and  imaging  display 
are  done  in  about  5  seconds,  which  is  fast  enough  for  clinical  experiments.  To  ensure  good  SNR, 
all  the  electronics  are  done  using  printed  circuit  boards. 
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Figure  7  shows  the  picture  of  the  entire  system  and  the  combined  probe.  Both  NIR  system  (top) 
and  ultrasound  system  (bottom)  are  mounted  on  a  hospital  cart.  The  combined  probe,  which 
houses  the  ultrasound  array  and  the  NIR  source  detector  fibers,  is  designed  to  be  hand  held  to 
scan  patients. 

C.  Phantoms 

0.5-0.6%  Intralipid  solutions  are  used  to  mimic  normal  human  breast  tissues  in  all  experiments 
and  the  corresponding  reduced  scattering  coefficient  fis  ranges  from  5  cm1  to  6  cm'.  The 

Intralipid  is  contained  in  a  large  fish  tank  to  set  up  approximately  a  semi-infinite  homogeneous 
phantom.  Small  semi-spherical  balls  (1  cm  in  diameter),  made  of  acrylamide  gel  ,  are  inserted 
into  Intralipid  to  emulate  lesions  embedded  in  a  breast.  The  reduced  scattering  coefficients  of 
the  gel  phantoms  are  similar  to  that  of  the  background  medium  {/Js  ~6  cm  ),  and  the  absorption 
coefficients  are  changed  to  different  values  by  adding  different  concentrations  of  India  ink  to 
emulate  high  contrast  (//a  =0.25  cm"1)  and  low  contrast  ( jia  =0.1  cm"1)  lesions.  Ultrasound 
scattering  particles  of  200  Jim  in  diameter  are  added  to  the  gel  phantom  before  the  gel  has  been 
formed. 

4.  Experimental  Results 

A.  Effects  of  missing  optical  sources  in  the  combined  probe 

A  series  of  experiments  was  conducted  to  estimate  the  optimal  hole  size.  Three  probe 
configurations  were  investigated,  which  were  (a)  no-hole,  (b)  2x2  cm2  central  hole,  and  (c)  3 
x3  cm2  hole  probes.  The  no-hole  probe  was  essentially  the  same  as  case  (b)  except  that  an 
additional  light  source  was  added  in  the  middle.  Figure  8  shows  reconstructed  NIR  images  for 
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on-center  targets  of  high  (jla  =  0.25  cm'1,  left  column)  and  low  contrast  (jua=0A  cm'1,  right 
column)  located  2.5  cm  deep  inside  the  Intralipid.  The  fitted  background  fda  and  jus  are  0.015 
cm'1  and  5.36  cm'1,  respectively.  With  the  target  depth  provided  by  ultrasound,  we  performed 
reconstruction  in  the  target  layer.  The  centers  of  the  voxels  in  this  layer  were  (x,  y,  2.5  cm), 
where  x  and  y  were  discrete  spatial  x-y  coordinates,  and  the  thickness  of  the  layer  was  1  cm. 
For  the  high  contrast  target  case,  there  are  no  important  differences  in  image  quality  associated 
with  different  probes  (Figure  8(a),  (c))  except  that  with  a  3  cm  by  3  cm  hole.  The  first  row  of 
Table  1  provides  measured  maximum  fJ.a  values  from  the  corresponding  images.  Because  of  the 
low  spatial  resolution  of  diffusive  imaging,  the  boundaries  of  targets  are  not  well  defined.  The 
maximum  value  is  a  better  estimation  of  reconstructed  target  fia .  From  no-hole  to  2  x  2  cm2,  the 

reconstructed  maximum  fia  decreases  slowly.  But  for  3x3  cm2,  the  maximum  /ia  drops 

suddenly  to  0.104  cm'1,  which  is  less  than  half  of  the  original  value.  Another  imaging  parameter 
we  measured  is  the  full  width  at  half  maximum  (FWHM)  of  the  corresponding  images.  Since  the 
image  lobes  were  elliptical  in  general,  we  measured  the  widths  of  longer  and  shorter  axes  and 
used  the  geometric  mean  to  estimate  FWHM.  The  results  are  shown  in  the  second  row  of  the 
Table  1,  and  the  FWHM  almost  increases  with  the  hole  size.  We  also  measured  image  artifact 
level,  which  was  defined  as  the  ratio  of  the  peak  artifact  to  the  maximum  strength  of  the  image 
lobe  and  was  given  in  decibels.  The  results  are  shown  in  the  third  row.  No  artifacts  were 
observed  in  the  images  of  no-hole  and  2x2  cm2  hole  probes.  However,  the  peak  artifact  level  at 
the  -14.3  dB  level  was  measured  in  the  image  of  the  3x3  cm2  hole  probe.  When  the  contrast 
was  low,  the  reconstructed  maximum  absorption  coefficients  and  measured  FWHMs  were 
essentially  the  same  for  the  no-hole  and  2x2  cm2  hole  probes.  However,  the  reconstructed 
maximum  value  dropped  to  60%  of  the  true  value  for  the  3x3  cm  probe.  The  artifact  levels 
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measured  in  images  of  three  probe  configurations  were  similar  and  were  worse  than  the  high 
contrast  case.  The  image  artifacts  are  related  to  the  reconstruction  algorithm.  When  the  target 
contrast  is  weak  and/or  the  SNR  is  low,  the  inversion  algorithm  produces  artifacts  around  the 
edges  of  the  images. 

For  shallow  targets  (here  we  set  the  target  depth  to  be  1.5  cm)  the  NIR  system  has  relatively 
poorer  performance.  This  is  due  to  less  source-detector  pairs  experiencing  the  existence  of  a 
shallow  absorber.  As  shown  in  Figure  9,  image  artifacts  are  obviously  worse  compared  with 

Figure  8.  However,  the  conclusion  about  the  hole  size  of  the  probe  remains  true.  Table  2  lists  all 

2 

the  measured  imaging  parameters  obtained  from  three  probe  configurations.  While  a  3  x  3  cm 
hole  is  somewhat  too  big  to  get  good  enough  results,  the  optimal  hole  size  is  in  the  neighborhood 
of  2x2  cm2 .  This  optimal  size  is  about  the  size  of  commercial  ultrasound  transducers. 

In  the  above  studies,  we  have  used  the  iteration  number  obtained  from  the  no-hole  configuration 
for  the  rest  of  the  configurations.  Ideally,  the  iteration  should  stop  when  the  object  function  (see 
Eq.  (2))  or  the  error  performance  surface  reaches  the  noise  floor.  However,  system  noise, 
particularly  coherent  noise,  was  difficult  to  estimate  from  experimental  data.  In  general,  we  have 
found  that  the  reconstructed  values  were  closer  to  true  values  when  the  object  function  reached 
about  5%  to  15%  of  the  initial  value  (total  energy  in  the  measurements).  Therefore,  we  have  used 
this  criterion  (-10%  of  the  initial  value)  for  the  no-hole  configuration.  Since  the  SNR  of  the  date 
decreased  with  the  increase  in  hole-size,  we  could  not  find  consistent  criterion  for  both  no-hole 
and  hole  data.  Therefore,  we  have  used  the  same  iteration  number  obtained  from  the  no-hole 
case  for  the  hole  configurations  and  the  comparison  was  based  on  the  same  iteration  number. 
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B.  Ultrasound-guided  NIR  imaging 

Three-dimesional  ultrasound  images  can  provide  3-D  distributions  of  targets.  With  the  a  priori 
target  depth  information,  the  optical  reconstruction  can  be  improved  significantly.  An  example 
is  given  in  Figure  10.  The  target  again  was  1  cm  diameter  gel  ball  of  low  (  //fl  =  0.1  cm'1)  optical 
contrast  and  was  embeded  at  approximately  (  0,  0,  2.5  cm)  inside  the  Intralipid  medium.  The 
fitted  background //fl  and  //’are  0.02  cm'1  and  5.08  cm'1,  respectively.  The  combined  probe 

shown  in  Figure  2  was  used  to  obtain  the  ultrasound  and  NIR  data  simultaneously.  Figure  10  (a) 
shows  an  A-scan  line  of  a  returned  ultrasound  echo  signal  received  by  one  ultrasound  transducer 
located  on  top  of  the  target.  As  acoustic  scatters  were  uniformly  distributed  in  the  target,  signals 
were  reflected  from  inside  the  target  as  well  as  from  the  surfaces.  The  reflected  signals  from  the 
front  and  back  surfaces  of  the  gel  ball  can  be  clearly  identified  in  the  echo  signal.  Based  on  the 
target  depth,  we  reconstructed  the  optical  absorption  coefficient  at  the  target  depth  only  (  1  cm  in 
thickness)  by  setting  the  purturbations  from  the  other  depths  equal  to  zero.  We  also  performed 
3-D  optical  only  reconstruction.  Figure  10  (b)  shows  the  recontructed  absorption  image  from  a 
3-D  optical  only  reconstruction  (layer  three  of  voxel  coordinates  (x,y,  2.5cm)  and  1  cm  thicknss), 
while  Figure  10  (c)  shows  the  reconstructed  image  of  the  corresponding  target  from  ultrasound 
guided  reconstruction.  For  optical  only  reconstruction,  the  algorithm  did  not  converge  to  a 
localized  spatial  region  and  the  image  contrast  was  poor.  The  measured  maximum  absorption 
coefficient  was  0.088  cm1,  which  was  close  to  the  true  value.  However,  the  measured  spatial 
location  of  the  maximum  value  was  (-1.6  ,  -1.2  cm)  which  was  too  far  from  the  true  target 
location.  With  the  a  priori  target  depth,  the  reconstruction  performed  at  the  target  layer  can 
localize  the  target  to  the  correct  spatial  position.  The  measured  maximum  absorption  coefficient 
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was  0.12  cm-1  and  its  location  was  (0,  0.4  cm)  which  was  very  close  to  the  true  target  location. 
This  example  demonstrates  that  a  priori  target  depth  can  significantly  improve  the  reconstruction 
accuracy  and  target  localization. 

In  addition  to  the  use  of  a  priori  target  depth  information,  we  can  also  use  target  spatial 
distribution  provided  by  ultrasound  to  guide  the  reconstruction.  We  performed  a  set  of 
experiments  with  two  targets  located  at  2.5  cm  depth  inside  the  Intralipid.  Each  target  is  a  1  cm' 
gel  cube  containing  ultrasound  scatters.  For  optical  properties,  they  both  could  be  high  contrast 
(f ua  =0.25  cm'1)  or  low  contrast  (jua  =0.1  cm'1),  but  had  the  same  reduced  scattering  coeficient 

as  the  background.  The  fitted  background  jua  and  ju's  are  0.017  cm'1  and  4.90  cm'1,  respectively. 

One  target  was  appproximately  centered  at  (-1.0,  -1.0,  2.5  cm),  while  the  other  was  at  (1.0,  1.0, 
2.5  cm).  The  distance  between  the  centers  of  the  two  targets  was  2.8  cm. 

Figure  11(a)  is  the  ultrasound  image  of  two  high  contrast  targets.  As  the  field  of  view  of  the 
ultrasound  system  was  nearly  a  3  cm  by  3  cm  square,  these  two  targets  appeared  at  diagonal 
comers.  The  measured  peak  positions  of  the  two  targets  were  (-0.6,  -1.0  cm)  and  (1.0,  1.0  cm), 
which  differed  from  the  true  target  locations  by  only  one  voxel.  The  low  contrast  of  the 
ultrasound  image  is  related  to  the  speckle  noise.  Since  our  ultrasound  array  is  sparse,  the 
imaging  quality  is  not  state-of-the-art  (see  more  discussion  in  Section  5).  The  NIR  image  of 
these  targets  was  obtained  simultaneously  and  is  shown  in  Figure  11(b).  The  reconstruction  was 
performed  at  the  target  layer  by  taking  advantage  of  target  depth  information.  A  total  of  123 
iterations  was  used  to  obtain  Figure  1 1  (b).  The  measured  peak  positions  of  the  two  targets  were 
(-1.4,  -1.0  cm)  and  (0.6,  0.6  cm),  which  were  one  voxel  off  from  the  true  target  locations  (-1.0, 
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-1.0  cm)  and  (1.0,  1.0  cm),  respectively.  The  corresponding  reconstructed  absorption 
coefficients  were  0.242  cm'1  and  0.251  cm'1,  which  were  close  to  true  values.  However,  the  two 
targets  were  almost  connected  to  each  other  and  their  spatial  localization  was  poor.  For  low 
contrast  targets,  the  ultrasound  image  is  shown  in  Figure  1 1  (c)  and  the  measured  peak  locations 
of  the  two  targets  were  (-1.0,  -0.6  cm)  and  (0.6,  1.0  cm),  which  differed  from  true  target 
locations  by  only  one  voxel.  The  corresponding  NIR  image  is  shown  in  Figure  1 1  (d)  and  the 
measured  peak  locations  of  the  two  targets  were  (-2.2,  -1.0  cm)  and  (0.6,  1.0  cm).  The  left  target 
was  off  the  true  location  by  three  voxels.  The  corresponding  reconstructed  absorption 
coefficients  were  0.063  cm'1  and  0.1004  cm'1  at  87  iteration  steps.  As  one  can  see,  the  target 
shape  and  localization  were  poorer  than  those  in  the  high  contrast  case.  In  addition,  an  artifact 
appeared  at  the  edge  of  the  image. 

From  the  co-registered  ultrasound  images,  we  obtained  spatial  distributions  of  the  two  targets 
and  specified  target  regions.  Figure  12  (a)  and  (c)  show  the  —6  dB  contour  plots  of  Figure  11  (a) 
and  (c).  Applying  the  same  reconstruction  scheme  to  these  specific  regions,  we  obtained  Figure 
12  (b)  and  (c)  in  1  iteration.  The  reconstructed  absorption  coefficients  were  0.2357  cm-1  and 
0.219  cm'1  for  the  two  high  contrast  target  cases,  and  0.123  cm'1  and  0.131  cm'1  for  the  low 
contrast  case.  We  can  see  much  better  improvement  in  the  low  contrast  target  case,  when 
comparing  Figure  12  (d)  with  Figure  11  (d).  This  example  has  demonstrated  that  target 
morphology  information  provided  by  ultrasound  can  be  used  to  guide  the  optical  reconstruction 
in  the  specified  regions.  ^ 

The  result  regarding  the  iteration  step  is  significant.  As  we  discussed  earlier,  there  is  no  known 
stopping  criterion  to  terminate  the  iteration  because  it  is  very  difficult  to  estimate  the  noise  level 
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in  the  measurements.  With  the  a  priori  target  depth  and  spatial  distribution  provided  by  co¬ 
registered  ultrasound,  we  can  obtain  an  accurate  optical  absorption  coefficient  in  one  iteration. 
Therefore,  no  stopping  criterion  is  needed  for  the  inversion  algorithms.  However,  this  result  will 
need  to  be  further  evaluated  with  more  samples  of  different  contrasts. 

5.  Discussion 

Commercial  ultrasound  scanners  use  1-D  probes  which  provide  2-D  images  of  x-z  views  of  the 
targets,  where  x  and  z  are  spatial  and  propagation  dimensions,  respectively.  Such  x-z  images 
cannot  co-register  with  NIR  images,  which  are  obtained  from  x-y  views  of  the  targets.  Our 
current  2-D  ultrasound  array  is  capable  of  providing  x-y  views  of  the  targets,  which  can  be  used 
to  co-register  with  NIR  images.  However,  the  array  is  sparse  and,  therefore,  the  image  resolution 
is  not  state-of-the-art.  Nevertheless,  its  spatial  resolution  is  comparable  to  NIR  imaging  and  can 
be  used  to  guide  NIR  image  reconstruction.  With  3-D  ultrasound  guidance,  only  one  iteration  is 
needed  to  obtain  accurate  absorption  coefficients.  This  result  is  significant  because  no  stopping 
criterion  is  necessary.  More  studies  with  a  variety  of  target  contrasts  and  locations  will  be 
performed  to  verify  this  result. 

We  have  purchased  a  2-D  state-of-the-art  ultrasound  array  of  1280  transducer  elements  and  we 
are  building  a  multiplexing  unit  for  our  64-channel  electronics.  In  addition,  the  new  2-D 
transducer  size  is  about  2x3  cm2,  which  is  in  the  neighborhood  of  the  optimal  hole  size  we  have 
found  through  this  study.  With  the  new  2-D  ultrasound  transducer,  we  will  be  able  to  obtain 
high-resolution  ultrasound  images  and  delineate  the  target  boundaries  with  finer  details  for 
optical  reconstruction. 
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Ultrasound  contrast  depends  on  lesion  acoustic  properties  and  NIR  optical  contrast  is  related  to 
lesion  optical  properties.  Both  contrasts  exist  in  tumors  but  the  sensitivities  of  these  two 
modalities  may  be  different.  It  is  possible  that  some  early  stage  cancers  have  NIR  contrast  but 
are  not  detectable  by  ultrasound.  It  is  also  possible  that  some  lesions  have  acoustic  contrast  but 
no  NIR  contrast  or  low  NIR  contrast.  However,  ultrasound  is  routinely  used  as  an  adjunct  tool  to 
X-ray  mammography,  and  the  detection  sensitivity  of  ultrasound  together  with  palpation  and  X- 

f . \ 

ray  mammography  N  98.4y  according  to  a  clinical  trial  reported  in  Ref  [1].  For  those  lesions 

that  are  seen  in  X-ray  and  ultrasound  images,  we  will  be  able  to  use  ultrasound  guidance  as 

reported  in  this  paper  to  quantify  the  light  absorption  as  well  as  other  optical  parameters  more 

accurately.  The  role  of  optical  imaging  is  to  add  more  specificity  to  the  ultrasonically  detected 
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lesions.  For  those  lesions  that  are  seen  by  X-ray  and  not  shown  in  ultrasound,  we  will  need  to 

A 

carefully  mark  the  corresponding  lesion  regions  in  ultrasound  images  and  reconstruct  NIR 
images  of  these  regions  as  well.  Optical  imaging  will  again  add  more  specificity  to  the 
conventional  imaging  method.  Since  X-ray  provides  projection  slices  with  no  depth  information, 
direct  X-ray  guidance  to  optical  imaging  may  not  yield  meaningful  results.  It  would  be  desirable 
if  we  could  obtain  sensitivity  of  optical  imaging  alone.  However,  light  scattering  is  a  main 
problem  that  prevents  the  accurate  and  reliable  localization  of  lesions. 

In  the  reported  phantom  studies,  we  assigned  zero  perturbations  to  the  regions  where  no  targets 
were  present.  In  clinical  studies,  we  plan  to  segment  the  ultrasound  images  and  specify  different 
tissue  types  as  well  as  suspicious  regions  in  the  segmented  images.  We  then  reduce  the 
reconstructed  optical  unknowns  by  assigning  unknown  optical  properties  to  different  tissue  types 
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as  well  as  to  suspicious  regions.  Finally,  we  reconstruct  the  reduced  sets  of  unknown  optical 
properties.  We  expect  more  accurate  estimation  of  reconstructed  optical  properties  and  fast 
convergence  speed,  as  reported  in  the  paper.  However,  it  is  still  too  early  to  judge  the  clinical 
performance  of  the  combined  method;  further  clinical  studies  are  needed. 

Probing  regions  of  the  banana-shaped  diffusive  photons  depend  on  source-detector  separations 
and  measurement  geometry.  For  a  semi-infinite  geometry,  the  probing  regions  extend  further 
into  the  medium  when  source-detector  separation  increases.  This  is  why  we  have  multiple 
source-detector  pairs  of  various  separations  to  detect  targets  at  variable  depths  from  0.5  cm  to  4 
cm.  Of  course,  it  is  hard  to  achieve  uniformity  sensitivity  in  the  entire  region  of  interest.  For 
example,  a  superficial  target  (~lcm  deep)  would  cause  very  strong  perturbations  when  it  is  close 
to  a  source  or  a  detector,  but  will  result  in  much  weaker  signals  when  it  is  located  deeper. 
Normalization  of  scattering  photon  density  waves  with  respect  to  the  incident  waves  makes  it 
possible  for  reconstruction  algorithms  to  handle  the  huge  dynamic  range  of  signals  and  to  detect 
a  target  as  deep  as  4  cm.  This  normalization  procedure  was  applied  to  the  reconstruction 
algorithm  used  to  obtain  the  reported  images. 

In  this  study,  the  target  absorption  coefficient  was  reconstructed  from  the  measurements.  Since 
the  target  jus  was  similar  to  the  background  fis ,  the  coupling  between  fia  and  jUs  in  our 

measurements  was  negligible.  We  also  did  experiments  with  gel  phantom  made  with  Intralipid 
of  the  similar  concentration  to  the  background  and  did  not  observe  perturbation  beyond  the  noise 
level.  Similar  reconstruction  studies  can  be  done  for  scattering  coefficients  as  well. 
Simultaneous  reconstruction  of  both  absorption  and  scattering  coefficients  is  also  possible. 
Since  the  eigenvalues  of  the  absorption  and  scattering  weight  matrixes  are  significantly  different. 
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good  regulation  schemes  are  needed  for  simultaneous  reconstruction.  This  subject  is  one  of  our 
topics  for  further  study. 

5.  Summary 

We  have  constructed  a  near  real  time  imager  that  can  provide  co-registered  ultrasound  and  NIR 
images  simultaneously.  This  new  technique  is  designed  for  improving  the  specificity  of  breast 
cancer  diagnosis.  Since  the  ultrasound  transducer  needs  to  occupy  the  central  region  of  the 
combined  probe,  a  series  of  experiments  has  been  conducted  to  investigate  the  effects  of  missing 
optical  sensors  in  the  middle  of  the  combined  probe  upon  the  NIR  image  quality.  Our  results 
have  shown  that  as  long  as  the  central  ultrasound  transducer  area  is  in  the  neighborhood  of  2  x  2 
cm2,  essentially  similar  reconstruction  results  as  those  of  no  missing  optical  sensors  in  the  middle 
of  the  combined  probe  can  be  obtained.  This  2x2  cm2  dimension  is  about  the  size  of  most 
commercial  ultrasound  phased-array  transducers.  When  the  central  missing  optical  sensor  area  is 
increased  to  3x3  cm2,  however,  the  reconstructed  values  are  obviously  lower  than  real  values.  If 
we  increase  the  iteration  steps,  artifacts  in  the  reconstructed  images  would  soon  become 
dominant. 

With  the  target  3-D  distribution  provided  by  co-registered  ultrasound,  significant  improvements 
in  algorithm  convergence  and  reconstruction  speed  have  been  achieved.  In  general,  a  priori 
target  depth  information  guides  the  inversion  algorithm  to  reconstruct  the  heterogeneities  at  the 
correct  spatial  locations  and  improves  the  reconstruction  speed  by  an  order  of  magnitude.  In 
addition,  a  priori  target  spatial  distribution  can  further  reduce  the  iteration  to  one  step  and  also 
obtain  accurate  optical  absorption  coefficients.  Given  the  fact  that  no  known  stopping  criterion 
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is  available,  this  result  is  significant  because  no  iteration  is  needed.  However,  this  result  will 
need  to  be  evaluated  with  more  samples  of  different  contrasts. 
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Table  1.  Parameters  of  reconstructed  images  for  deep  high  (fia  =0.25  cm  )  and  low 


(fj.a  =  0.1  cm’1)  contrast  targets.  The  fitted  background  fia  and  jus  are  0.015  cm  1  and  5.36  cm 

A 

respectively.  //a(max)  is  the  measured  maximum  value  of  the  reconstructed  absorption 

coefficient  map,  FWHM  is  defined  as  the  geometric  mean  of  the  widths  measured  at  longer  and 
shorter  axes  of  the  elliptical  image  lobe. 


Probe  type 

No-hole 

2  x2 

3x3 

High  contrast 

A 

Vain ax)  (cm'') 

0.251 

0.234 

0.104 

FWHM  (cm) 

1.85 

1.91 

2.44 

Artifacts  (dB) 

Background 

(-22dB) 

background 

-14.3 

Low  contrast 

A 

V  aimax)  (Cm‘‘) 

0.105 

0.111 

0.064 

FWHM  (cm) 

2.02 

1.83 

2.16 

Artifacts  (dB) 

-6.90 

-8.10 

-5.65 
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Table  2.  Parameters  of  reconstructed  images  for  shallow  high  (jua  =0.25  cm"1)  and  low 
(jua  -  0.1  cm'1)  contrast  targets.  The  fitted  background  fia  and  ju’s  are  0.015cm'1  and  5.35  cm1, 


respectively. 


Probe  type 

No-hole 

2  x2 

3x3 

High  contrast 

A 

/^u(max)  (Cm') 

0.250 

0.194 

0.118 

FWHM  (cm) 

1.32 

1.61 

2.08 

Artifacts  (dB) 

-7.98 

-12.7 

-9.76 

Low  contrast 

A 

Main**)  (Cm1) 

0.100 

0.091 

0.042 

FWHM  (cm) 

1.88 

2.11 

3.17 

Artifacts  (dB) 

-6.25 

-7.44 

-0.65 
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Figure  Captions 

Figure  1.  Schematic  arrangement  of  NIR  source  and  detector  fibers  on  the  probe.  Small 
solid  circles  are  source  fibers  and  bigger  solid  cycles  are  detector  fibers. 

Figure  2.  Picture  of  an  experimental  probe.  An  ultrasound  array  of  8  x  8  =  64  transducers 
occupies  the  central  3x3  cm2  area,  and  12  dual  wavelength  source  fibers  and  8  detector 
fibers  are  deployed  at  the  periphery. 

Figure  3.  Schematic  of  the  NIR  frequency  domain  imaging  system.  The  modulation  frequency 
is  140MHz.  The  12  dual-wavelength  source  channels  are  switched  on  sequentially  by  a  PC  and 
8  detector  channels  receive  signals  in  parallel. 

Figure  4.  (a)  Log(plfl  Aa/} )  vs.  distancev  after  calibration,  (b)  phase  (j)ap  vs.  distance  pa/) 

after  calibration. 

Figure  5.  Schematic  of  our  ultrasound  scanner.  64  ultrasound  transducers  are  connected  to  64 
parallel  transmission  and  reception  channels.  The  transmission  part  consists  of  64  high  voltage 
pulsers,  which  can  be  controlled  by  computer  generated  delay  profiles.  The  reception  part 
consists  of  64  two-stage  amplifiers  and  A/D  converters. 

Figure  6.  Ultrasound  subarray  scanning  configuration. 

Figure  7.  Picture  of  our  combined  system.  NIR  system  (top)  and  ultrasound  system  (bottom)  are 
mounted  on  a  hospital  cart. 
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Figure  8.  Reconstructed  NIR  images  of  deeper  targets  (2.5  cm  in  depth,  1  cm  in  diameter,  and 
fitted  background  //u  and  //’are  0.015  cm'1  and  5.36  cm’1,  respectively).  The  left  column 
corresponds  to  images  of  a  high  contrast  target  (//a  =  0.25  cm'1)  obtained  from  different  probe 
configurations,  while  the  right  column  corresponds  to  images  of  a  low  contrast  target 
( na  =0.1  cm'1).  Each  row  is  related  to  a  specific  hole  size:  1)  No  hole  (a  and  b);  2)  2  x2cm  2  (c 

and  d);  and  3)  3  x  3  cm2  (e  and  f). 

Figure  9.  Reconstructed  NIR  images  for  shallow  targets  (1.5  cm  in  depth,  1  cm  in  diameter,  and 
fitted  background  jua  and  //’are  0.015  cm'1  and  5.36  cm"1,  respectively).  The  left  column 
corresponds  to  images  of  a  high  contrast  target  ( jua  =0.25  cm'1),  while  the  right  column 
corresponds  to  images  of  a  low  contrast  target  (//0  =0.1  cm'1).  Each  row  is  related  to  a  specific 
hole  size:  1)  No  hole  (a  and  b);  2)  2x2  cm2  (c  and  d);  and  3)  3x3  cm2  (e  and  f). 

Figure  10.  A  deep  target  (2.5  cm  in  depth,  1  cm  in  diameter)  of  low  optical  contrast 
(jua  =0.10 cm'1,  and  fitted  background  //fl  and  //’are  0.02  cm'1  and  5.08  cm'1,  respectively),  (a) 

A-scan  line  of  the  reflected  ultrasound  pulse-echo  signal  indicating  the  target  depth,  (b) 
absorption  image  of  the  low  contrast  target  obtained  from  optical  only  reconstruction,  (c) 
ultrasound  guided  reconstruction  at  target  depth. 


Figure  1 1 .  Simultaneously  obtained  ultrasound  and  NIR  absorption  images.  The  fitted 
background  fla  and  //’are  0.017  cm"1  and  4.90  cm'1,  respectively,  (a)  ultrasound  and  (b)  NIR 


absorption  image  of  two  high  contrast  targets  (target  //a  =0.25  cm'1),  (c)  ultrasound  and  (d) 

NIR  image  of  two  low  contrast  targets  (target  jiu  =  0.10  cm'1).  In  both  high  and  low  contrast 
cases,  the  two  targets  were  located  at  2.5  cm  depth. 

Figure  12.  (a)  and  (c)  are  -6  dB  contour  plots  of  ultrasound  images  shown  in  Fig.  11  (a)  and  (c). 
The  outer  contour  is  -6  dB  from  the  peak  and  the  contour  spacing  is  1  dB.  (b)  and  (d)  are 
corresponding  NIR  absorption  maps  reconstructed  in  target  regions  specified  by  ultrasound. 
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